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THE POISSON PASCAL DISTRIBUTION’ 


S. K. Karti? anp JoHn GuRLAND? 


Towa State University of Science and Technology 
Ames, Iowa, U: S. A. 


1. INTRODUCTION 


Elementary distributions such as the Poisson, the Logarithmic 
and the Binomial which can be formulated on the basis of simple 
models have been found to be inadequate to describe the situations 
which occur in a number of phenomena. The Neyman Type A (ef. 
Evans [5]), the Negative Binomial (cf. Bliss and Fisher [3]), and the 
Poisson Binomial (cf. McGuire et al. [8]), which combine two of the 
elementary distributions through the processes of compounding and 
generalizing (cf. Gurland [7]), have been fitted with varying degrees of 
success to data from a number of biological populations. The aim of 
this paper is to study what may be called the Poisson Pascal distribu- 
tion which includes the Neyman Type A and Negative Binomial as 


particular limiting cases and serves as a natural complement of the 
Poisson Binomial. 


2. FORMULATION OF POISSON PASCAL DISTRIBUTION: 


Consider the well-known example of the egg masses and the larvae 
illustrating the process of generalizing (cf. Neyman [9]). Suppose that 
the egg masses in a field have a Poisson distribution with probability 
generating function (p.g.f) exp. {A(z — 1)} and that the survivors 
within an egg mass have a Pascal‘ distribution with p.g.f. (¢ — pz)~*, 
p>0,q=1+p,k> 0. For the sake of references, we note that, 
if g(z) is the p.g.f. of a distribution, then the coefficient of z* in the 
Taylor's expansion of g(z) yields the probability of count z. Then, 
by following the arguments in Gurland [7], we obtain for the distribu- 


tion of the survivors in the field, the Poisson Pascal distribution 
with p.g.f. 


1This research was sponsored in part by the Air Force Office of Scientific Research. 

2Current address: The Florida State University, Tallahassee, Florida. 

3On leave at Mathematics Research Center, U. S. Army, University of Wisconsin. 

‘Feller [6] refers to the distribution with p.g.f. (¢ — pz)“*, p > =1+7,k =1,2,3°°* ete. 


However, we will use the work Pascal as synonymous with Negative Binomial and use it in a place 
of the latter since it is shorter. 
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exp — pz) *— 1]}, K>0, p>O0, g=1t+p. (I) 


It is believed that most of the deviation of the distribution of the 
survivors from the simple distributions mentioned above stems from (i) 
the complex structure of survivors within an egg mass and (ii) the 
movement of survivors from place to place. Since the Negative Bi- 
nomial distribution has been found to be very useful in fitting data 
involving this type of heterogeneity, it is reasonable to suppose that 
the Poisson Pascal will give a better description of the population of the 
survivors in a field. 

The Poisson Pascal distribution can also be looked upon as the 
result of compounding a Pascal distribution with p.g.f. (q — pz)" by 
taking k, to behave as kx where k is a positive constant and x a Poisson 
random variable with p.g.f. exp {A(z — 1)}. 


3. SOME PROPERTIES OF THE POISSON PASCAL DISTRIBUTION 


The limiting forms that the Poisson Pascal distribution takes as the 
parameters take on extreme values are given in Table 1. It is to be 


TABLE 1 


Some Limitine Forms oF THE Potsson Pascal DISTRIBUTION 
No. Limits Taken Name and p.g.f. of the limit 
1 k- o,p—0 Neyman Type A, exp {Afexp (Ai(z — 1) — 1]} 
pk =i 
2 Negative Binomial, (¢ — 
= ky 
3 Poisson, exp {\i(z — 1)} 
kp =X 


noted that the Neyman Type A and the Pascal distributions are among 
the limiting forms. 

The flexibility of the Poisson Pascal was compared quantitatively 
with that of the Neyman Types A, B, C, Pascal, and Poisson Binomial 
by evaluating the relative skewness and kurtosis of each for fixed 
mean k,p, and variance k,p,(1 + p,) using the indices x,3;/k,p? and 
kta)/kip} respectively along the lines of Anscombe [1]. The range of 
numerical values of these indices for each of the foregoing distributions 
is shown in Table 2. It is apparent that the Poisson Pascal covers the 
entire range of distributions from Neyman Type A to Pascal with 
respect to skewness and kurtosis. Also, the ranges of these ratios for 


q 
eg 
‘ 
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TABLE 2 
CoMPARISON OF SKEWNESS AND Kurrtosis oF CERTAIN DISTRIBUTIONS 
No. Name Range of Skewness Range of Kurtosis 

1 Neyman Type A 1 1 
2 Neyman Type B 9/8 27/20 
3 Neyman Type C 6/5 8/5 - 
4 Pascal 2 6 
5 Poisson Pascal (1,2) . (1,6) 
6 Poisson Binomial (0,1) (0,1) 


the Poisson Pascal and the Poisson Binomial are disjoint. Since their 
p.g.f.s have the common form 


exp {A[(g — pz)“ — (2) 
we observe that the distribution with (2) for p.g.f. wherein A > 0, 
q = 1+ p,p > 0whenk > Oand —1 < p < 0 when k is a negative 
integer, covers a very wide range of distributions. As an aid in com- 
puting the values of the ratios for the samples to obtain an idea as to 
how close the sample distribution is to the various distributions men- 
tioned in Table 2, formulae to compute the factorial cummulants x, ;; 
using sample moments are given in Appendix A. 

Since the first two frequencies were large in the sets of data to 
which this (and similar) distributions were fitted, the ratio of the first 
two frequencies were compared for some of these distributions. It 
can be easily shown that the value of this ratio for the Poisson Pascal 
distribution lies between the ratios for the Neyman Type A and the 
Pascal distributions. The ratios are given for brevity in Table 3. 


TABLE 3 


CoMPARISON OF THE RATIO OF THE FirsT TWO FREQUENCIES WITH MEAN 
AND VARIANCE FIXED AS kiP; AND k,P,(1 + p:) 


No. Distribution Ratio of Frequencies 
1 Neyman Type A kip; exp (—7p:) 
2 Neyman Type B 2k, {1 — exp (—3p:)(1 + 3p:)} (pig) 
3 Neyman Type C 3ki{2p, exp (—4pi) + 2[—2 exp (—4pu) 
+ 64p? + /4p} 
4 Pascal + pr) 
5 Poisson Pascal kipi(l + pi)/(k + 
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4. FITTING POISSON PASCAL AND EFFICIENCY OF 
METHODS OF ESTIMATION 


Since obtaining maximum likelihood estimates is very cumbersome, 
ad hoc methods were used to estimate the parameters. When the mean 
and the variance were large, use was made of the method of the first 
three moments. When they were moderate and the proportion of the 
zero frequency large, use was made of the method of the first two 
moments and the proportion of the zero frequency. When the first 
two frequencies were large in comparison with the remaining, the 
method of the first two moments and the ratio of the first two fre- 
quencies was used in estimation. The equations for estimation are 
given in Appencix C. The fit of this distribution to the data of Beall 
and Rescia [2] are given in Tables 4 and 5. For the sake of reference, 
the fit of a generalization of the Neyman Type A as given by the authors 
of these data are also given alongside. The relatively good fit of the 
Poisson Pascal is apparent. 


For obtaining a comparison of the various methods of estimation, 
it was decided to compute the efficiency function 


TABLE 4 
Fit oF THE OBSERVED FREQUENCY OF LesPEDEZA CapiTaTa, TABLE V oF [2] 
Expected frequency 
Observed due to Poisson Pascal Expected frequency 
Plants Frequency (Method of Moments) as in [2] 

0 7178 7185.0 7217.6 
1 286 276.0 218.6 
2 93 94.5 105.5 
3 40 41.5 50.9 
4 24 20.2 24.5 
5 7 10.4 11.8 
6 5 5.6 5.7 
7 1 3.1 2.8 
8 2 1.3 
9 1 1.0 6 

10 2 6 3 

11+ 1 3 4 

9.58 42.97 

Degrees of 
Freedom _ 8 9 


ie 
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E = 1/(Generalized variance X Information determinant) 
(cf. Cramer [4], pp. 489-497). (3) 


If we denote the parameter vector (A, p, k) by (A; , Az , As) and the 


set of statistics used by (¢, , t2 , 3), we get the expression for the general- 
ized variance of the estimates as 


TABLE 5 


Fit or OBSERVED FREQUENCY OF LEPTINOTARSA Deanna, 
TaBLeE III oF [2] 
Expected frequency due to 
Poisson Pascal, (method of Expected 
two moments and first frequency 
Insects Observed frequency as in [2] 

0 33 33.0 39.5 

1 12 9.8 6.0 

2 5 7.4 4.9 

3 6 5.5 3.4 

4 5 4.0 3.2 

5 0 2.9 2.5 

6 2 2.1 2.0 

7 2 1.5 1.6 

8 2 ye | 1.3 

9 0 8 1.0 
10 8 
11+ 2 1.3 3.3 
= 6.88 13.75 

Degrees of 
Freedom —_ 8 9 
a(n; ) 2 
| V(t, ~ ds) (4) 


where 7; , 72 , 7s are the functions of ), , \, and A; estimated consistently 
by t,,4,%. A proof of this is given in Appendix B. Since evaluating 
the covariance matrix V(t, , t, , t;) and the derivatives of ¢, , t, and & 
follows from the regular statistical techniques, no elaboration need be 
made here. A formula for the information determinant (cf. Shenton 
[10}) is 


= 
cay 
4 


532 BIOMETRICS, DECEMBER 1961 


1 dP, oP, 1 dP, oP, 
Ai OA, dA, a, AA, OAs 


= 1 aP, aP, 1 aP,aP, 1 OP, oP, 
P. Ae OA, P, Or2 z Or2 OA3 (5) 


@P, aP, 1 oP, aP, 1 oP, aP, 


The principal problems in evaluating J for a value of (A, p, k) therefore 
are (i) to evaluate the various P, and the derivatives of P, and (it) 
to determine the number of terms to be used in summing the infinite 
series. To obtain P, , let 


g(2) = exp {Al(q — — 1)} (6) 


and 
= (q — pe)*. 
By differentiating (6) successively we get 


g'(@) = Ag@h'@), (7) 
and 


Set z = 0 in equations (6), (7) and (8) and observe that g“(0) = 2! P, 
and h‘*(0) = zx! x, where 


! 
= (9) 


is the probability of z in the Pascal distribution with A(z) as p.g.f. 
Then we have the recurrence formulae 


P, = exp. {r[(g)* — 1)} (10) 
and 


Pon = = {Ee +1 iP}, (11) 


i=0 


which can be repeatedly used to evaluate any P, . 


The various derivative are given by 


oP, _ oP, 


(x + Pi 
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1 
= — (x + (13) 
and 
oP, _ 1 [2 |} 
a; ok x! Lar 9” 
= log (q — 
(p/9)' = {ae — + — — 
1 
— (log 9){q(z + 1)P..1 — (14) 
for all z. 
TABLE 6 
EFFICIENCY OF THE METHOD OF THE First THREE MOMENTS FOR 
THE Poisson PascaL 
k 
r Pp 1 3 5 1.0 2.0 
A .84 .82 .82 -76 
1.0 .26 .25 .24 .22 
2.0 13 12 13 
5 .90 .82 77 .67 
3 .59 .58 49 .35 
5 .5 .46 44 Al 34 .22 
1.0 .26 25 18 
1.0 81 .83 .82 75 
1.0 3 .59 59 .55 .46 .28 
1.0 .46 .44 .40 ol 
1.0 1.0 27 .25 23 15 -05 
5.0 94* 99* .78* 58* 
5.0 62* 56* 11° .01* 
5.0 5 48* 38* .20* 04* .00*. 
5.0 29* 17° .06* 02* .00* 


: 
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To determine the number of terms we use the following rule: 

Let 7;;(n) denote the (n + 1)th term in the series involved in the 
(i, j)th term of matrix (5). Let S;;(n) = >o*., 7T.,(r). Compute 
S.;(n) and >>,; T?,(n)/S?,(n) for n = 1, 2, --- et cetera successively 
till a value of n is reached for which 


Tisn)/Si(n) < 10™*. (15) 


It is clear from (15) that | 7,;(n) |/| S,;(n) | <10~* for each 7 and j. 
If the series converge faster than a geometric series with common 
ratio less than 0.9 and this convergence starts before the value of n is 
reached for which (15) is satisfied, the calculated efficiency will be 
correct to three significant figures. If the significant figures do not 
cancel out, this should yield the efficiencies computed therefrom, correct 
to three decimal places. When the inequality (15) was not satisfied 
for values of n < 20, the partial sum of the first twenty terms was 
taken as the value of the series since evaluating the terms of the series 


TABLE 7 


EFFICIENCY OF THE METHOD OF THE First Two MoMENTs AND THE First 
FREQUENCY FOR THE Poisson PAscAL DISTRIBUTION AT A = 0.1 


k 
N p 1 3 5 1.0 2.0 
99 .98 98 .98 98 
3 94 95 .94 93 
5 .90 .90 91 91 89 
r 1.0 82 83 84 84 83 
2.0 74 76 78 84 78 
5 A 1.00 1.00 98 96 95 
5 3 93 94 93 .92 87 
5 5 91 91 .90 89 82 
5 1.0 83 85 85 82 73 
5 2.0 75 78 81 
1.0 A 1.00 99 97 96 94 
1.0 2 94 96 94 93 87 
1.0 5 92 92 91 98 81 
1.0 1.0 84 85 86 82 72 
1.0 2.0 88* 87* 87" .87* 83* 
5.0 .99* 89" 98" 
5.0 3 97* 1.00" 98* 
5.0 5 95* .97* 98" 99* 96" 
5.0 1.0 .90* 96" 97* 
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for n larger than 20 is very time consuming. The efficiency when n 
was restricted to 20 is marked with an asterisk to indicate that they 
are less likely to be correct to three decimal places. When the efficiency 
so computed was larger than one (due to the inaccuracy in computing 
the information determinant), the corresponding cell in the efficiency 
table is left blank. 

The efficiency of the method of moments is given for certain values 
of (A, p, k) in Table 6. The efficiency of the method of the first two 
moments and the first frequency is given in Table 7 and that of the 
method of the first two moments and the ratio of the first two fre- 
quencies in Table 8 for the same values of (A, p, k). 

It is apparent that the method of the first two moments and the 
ratio of the first two frequencies has high efficiency and is superior to 
the other two methods when A is small (and consequently the first 
two counts account for a large proportion of the observed frequencies). 
Also the method of the first two moments and the first frequency is 
highly efficient when A is moderately large. When X, p or k approaches 


‘TABLE 8 


EFFICIENCY OF THE METHOD OF THE First Two MOMENTS AND THE RATIO OF 
THE First TWO FREQUENCIES FOR THE Poisson Pasca aT A = 0,1 


k 
r Pp 1 3 5 1.0 2.0 
| ‘1.00 1.00 1.00 1.00 1.00 
3 .98 .99 .99 .99 .98 
5 .99 .98 .98 .97 94 
1.0 94 .93 91 
2.0 .87 .88 .89 91 -99 
5 a 1.00 .99 1.00* 1.00 .99 
5 3 .98 1.00 .99 -96 
5 1.0 99 .97 .96 -92 
5 2.0 94 .94 .93 .89* .80* 
1.0 .98 1.00 1.00 1.00 
1.0 3 .99 1.00 .99 .99 95 
1.0 5 1.00 .98 .96 91 
1.0 1.0 95 94 93 .88 .79 
1.0 2.0 .88* .89* 
5.0 a 99* .95* 91* 
5.0 3 1.00 99* 98* .89* 76" 
5.0 5 1.00 96* 91* .80* 65* 
5.0 1.0 .95* 88* 79° 84* 1.00* 
5.0 2.0 .90* 98* 


1 

| 
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infinity, it can be shown by calculus that the method of moments is 
much superior to the other two but since the efficiency of each of these 
methods tends to zero for such values, this has little significance. 


5. CONCLUSIONS 


On the basis of the properties discussed in Section 3 and the fitting 
in Section 4, we observe that the Poisson Pascal distribution acts as a 
bridge between the Neyman Type A and the Negative Binomial dis- 
tributions and may be used with advantage when the latter distribu- 
tions are inadequate to represent the population accurately. 

From the tables of efficiency, it is clear that in the region of tabula- 
tions, at least one of the ad hoc methods of estimation suggested above 
has high efficiency. It is believed that in practice, (A, p, k) will not 
be far beyond the region of tabulation and that one of these methods 
can be used without too much loss of information. Techniques for 
choosing one of the many ad hoc methods on the basis of the sample 
will be discussed in a future paper. 
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APPENDIX 


A. Formulae to compute factorial cumulants using moments about the origin: 


We first obtain formulae to compute the first four factorial moments 
His) ,7 = 1, -+- , 4 using moments about the origin uv; and then obtain 


x 
: 
: 
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formulae to evaluate the first four factorial cumulants «,;; i = 1, --+ , 4 
using these factorial moments. 

As for the first objective we note that u;;,;=E{z(r—1)---(e—i+1)}, 
i = 1,2,---. By expanding the product within the expectation sign 
and using the elementary properties of the expectation operator, we get 

Huy = = Mis) = — + Qui, 
and 
Mia) = — Gus + — (16) 
As for the latter, we observe that if u(t) and y(t) denote the factorial 
moment generating function and the factorial cumulant generating 
function, then y(t) = log u(t). On differentiating the equation 
successively with respect to ¢ at ¢ = 0 and noting that «,,,; = ¥‘(0), 
we have 
Ka) = Mi» Ki2) = — Bi 
Kis) = Mis) — + 
and 
Kia) = Bia) — — + = (17) 
B. To prove formula (4): 

Since 7; , 72 , 73 are functions of ), , A, and A; , let us write them more 
explicitly as , Az , As), , Az As) and 73 (Ay, Az , As) Fespectively. 
If \, , X42 , As are the estimates of A, , Az , As , then using the statistics 
t, , tg , ts which estimate the 7’s consistently, we have 

i; = » As); i= 1, 2, 3. 
Hence, we have the asymptotic relations 
0 0 Or; 
which can be rewritten as 


— Ag | T. 


The generalized variance G of i, , i, , X3 is then given by 
G = | EQ, — A, — Any As — As)", — Ad, Ae — Az, As — As) | 


= 1 | 


which is formula (4). 


1 
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C. Equations for estimation: 


We give below the equations for estimating the parameters for the 
various methods mentioned in Section 4. Their derivations are omitted 
for brevity. 


i. Equations for the estimation of parameters using the first three 
moments are: 


Risky — (Ria)? 
= + 1), (19) 
A = &/(kp). (20) 


ii. Equations for the estimation of parameters using the first two 
moments and the proportion P, of zeros are: 


p log + (i - + (ie r) log(1+p), (21) 
k= 421 = (22) 


iii. Equation for the estimation of the parameter p using the first 
two moments and the ratio P,/P, of the first two frequencies is 


k and ) are then estimated by using equations (22) and (23). 


> 
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SOME RANK SUM MULTIPLE COMPARISONS TESTS 


Rosert G. D. Stee. 


Institute of Statistics, North Carolina State College 
Raleigh, North Carolina, U.S. A. 


1. SUMMAR™ 


Rank sum tests for multiple comparisons and suitable to the 
completely random design with equal sample sizes are discussed. 
Significance tables and some facts about the joint distribution of rank 
sums are given. An example illustrating test procedures and making 
use of the tables is presented. 


2. INTRODUCTION 


A number of nonparametric tests based on ranks have been pro- 
posed for the comparison of treatments in a completely random design. 
For example, we have the Wilcoxon-Mann-Whitney test [21, 10], 
basically a two-sample test with a per-comparison error rate. 

Also, Kruskal and Wallis [9] have proposed a rank test which is an 
analogue of Snedecor’s F-test. This test provides evidence concerning 
the presence of real differences but is of limited use in locating them. 

Steel [16, 17] has presented rank tests for comparing treatments 
against control and for all pairwise comparisons. Both of these tests 
use experiment-wise error rates. | 

Pfanzagl [13], as part of a more general theory, has discussed a 
two-step nonparametric decision process based on ranks, for testing 
the null hypothesis that k samples come from the same population and, 
if this is rejected, for deciding which one of the samples comes from a 
different population. No tables are given but it is suggested that they 
might be obtained by random sampling. It is also shown that the 
limiting distribution of the multivariate criterion is multinormal. 

The per-comparison error rate test is sometimes criticized, par- 
ticularly when all possible paired comparisons are made, because it 
will almost certainly lead to false declarations of significance when the 
experiment includes many treatments and if customary significance 
levels are used. It is also deemed inappropriate when the experiment 
is considered to be the conceptual unit. 

The experimentwise error rate test is sometimes criticized because it 
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requires such a large difference to be declared significant that it becomes 
difficult to detect any but the largest real differences when customary 
significance levels are used. Also, it may be that the individual com- 
parison is considered to be the conceptual unit. 

A brief discussion of these error rates is given by Steel [18] in re- 
sponse to Biometrics Query 163. 

Choice of a definition of error rate in the conduct of a particular 
experiment seems somewhat less crucial when it is realized that we can 
compute the significance level for a particular definition of error rate 
from knowledge of the chosen significance level for any other definition 
of error rate. This is not generally a simple computation unless the 
comparisons are independent. In the case of p independent comparisons, 
if a’ and a@ are the experimentwise and per-comparison error rates 
respectively, we have the relation: 1 — a’ = (1 — a)’. 

When comparisons are not independent, computation of comparable 
significance levels for different definitions of error rate depends upon 
the extent of the dependence and the nature of the multiple comparisons 
test. No individual is likely to perform such a computation for a 
single experiment. Thus a table needs to be prepared for comparing 
the customary significance levels for differently defined error rates. 
For tests based on an underlying normal distribution, this has been done 
fairly extensively by Harter [4, 6]. 

The experimenter may try to meet the usual criticisms of per- 
comparison and experimentwise error rates by choice of a non-standard 
significance level or an alternative test procedure. Presently, tables 
of significant values for such levels do not appear to be available for 
experimentwise error rates; in the case of alternative tests, several are 
available, including the Newman-Keuls [11, 7] procedure and Duncan’s 
[1,5] new multiple range test, which are sequential in nature. 

This paper is concerned with rank tests, in particular with tables for 
an all-pairs-of-treatments test with an experimentwise error rate, and 
the use of these tables for a fixed rank sum test, an analogue of Tukey’s 
w-procedure (for example, see Steel and Torrie [19]), and for two multiple 
rank sum tests, analogues of the Newman-Keuls procedure and of 
Duncan’s test. Table 2 is used in the first two cases, Table 3 in the 
last case. 


3. CONSTRUCTION OF TABLES 


The proposed tests call for rank sums and their conjugates computed 
for all pairwise comparisons of treatments. The minimum of each 
rank sum and conjugate is used, the set of minima providing a multi- 
variate rank sum test criterion. These sums are compared with a single 
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tabuiaved value for the fixed rank sum test and with several values for 

the multiple rank sum tests. Table 2 provides critical values for the 

analogues of Tukey’s and the Newman-Keuls tests; Table 3 provides | 
for the analogue of Duncan’s test. 

Methods for constructing probability tables and limited tables have 
been presented earlier [16, 17]. Construction of exact tables of any 
extent was beyond the computing facilities available. However, some 
machine time was available and this was used for some sampling ex- 
periments. 

It was originally intended to ignore the discrete nature of the data 
and to use the Kolmogorov-Smirnov [8, 14] one-sample test to determine 
the sample size necessary to attain a certain precision in the constructed 
tables. However, available computing facilities limited sampling to 
values of k = 3 and 4 (2 and 3 treatments when one was control) and 
n = 4 (1) 10. In addition, samples were obtained for k = 5, n = 4, 
5 and k = 6,n = 5. The number of permutations obtained for each 
case was either 5000 or 6000. These tables were used only for checking 
purposes against the few exact distributions available, k = 3 and 
n = 3 (1) 6, and against approximations used in constructing these 
and earlier tables. 

It was assumed that the various multivariate rank sum criteria 
are distributed approximately as multinormal distributions having 
mean vectors and variance-covariance matrices as given in the appendix. 
(Fraser’s [3] vector form of the Wald-Wolfowitz-Hoeffding theorem 
does not apply since the || C,.(7, j) ||'s of Fraser do not exist for the 
test criteria used here.) _ 

On this assumption, one naturally proceeds to base computation 
on presently available tables. Tables for Tukey’s and Duncan's tests 
are the obvious choice for all-pairs tests. These tests are based on a 
multinormal distribution with p? = n,n,;/(n, + n;)(m, + n;) where the 
present distribution calls for p? = njn;/(n, + n; + 1)(m +n; + 1), 
a small difference. The appropriate tables are, then, tables of the 
Studentized range with known variance, that is, infinite degrees of 
freedom. Table 22 of Pearson and Hartley [12] is such a table, gives 
percentage points of .10, .05 and .01, is appropriate for the first two 
tests, and was used in computing Table 2. Corrected tables for Duncan’s 
test have been computed by Harter [5] and this table was used in 
eomputing Table 3, also for percentage points of .10, .05 and .01. 

Table 2 was constructed by taking the integral part of u — to/ V2, 
unless the decimal fraction was > .9, in which case the next higher 
integer was tabulated, where ¢ was obtained from the distribution of 
w/o, w = range, Table 22 of Pearson and Hartley [12]. Since rank 


{ 
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sums are essentially differences, it is necessary to introduce +/2 into 
the denominator as shown. Tabulated rank sum values for a = .10 
differed in only two cases from values obtained by sampling. In par- 
ticular, for k = 4, n = 4, no value is significant by sampling; for k = 3, 
n = 6; a rank sum of 25 is significant by sampling whereas 26 is not. 
Values for a = .05 differed in no case. Values for a = .01 ran lower 
than those obtained by sampling, the difference increasing with n to a 
value of two in three cases. Hence, it is reasonable to conclude that 
tabulated values of the rank sum are conservative (low) for a = .01. 

The first attempt to construct Table 3 led to values which tended to 
run high for a = .10, correct for a = .05, and low for a = .01, relative to 
values found by sampling. For this reason, tabulated values are of 
u — to/~/2 decreased by unity for a = .10, as computed for a = .05, 
and increased by unity fora = .01. On this basis, tabulated values 
for a = .10 appear to be low, hence conservative, when not in agreement 
with sampling results; in particular, 7 out of 19 tabulated values are 
one unit low. For a = .05, 3 out of 19 values are one unit high. For 
a = .01, 5 out of 19 values are one unit low, with three of these being 
forn = 5. 

Tables have already been constructed [16], using Dunnett’s [2] 
tables, for the treatments against control test. These tables agree well 
with the sampling results. In no case is there a difference of more 
than one in the value of the test criterion, with the tables most often 
giving the conservative (lower) value. 


TABLE 1 


Finat Weicuts or Cuicxs at Weexs (GRAMS) For VARIOUS 
Sources oF PROTEIN SUPPLEMENT 


H L Sb Sf M 
Horse- Linseed Soybean Sunflower Meat C 
bean Oil Meal Oil Meal Seed Oil Meal Meal Casein 


: 
| 
179 309 243 423 325 368 
160 229 230 340 257 390 
136 181 248 392 303 379 
: 227 141 327 339 315 260 
. 217 260 329 341 380 404 
‘ 168 203 250 226 153 318 
2 108 148 193 320 263 352 
: 124 169 271 295 242 359 
143 213 316 334 206 216 
‘3 140 257 267 322 344 222 


RANK SUM MULTIPLE COMPARISONS 543 


4. USE OF TABLES 


To illustrate the use of the tables, the data in Query 60 of Biometrics 
(15) are used. These are presented in Table 1. Since the test is presently 


unavailable for unequal sample sizes, only the first ten items in each 
treatment are used. 


The following set of minimum rank sums is obtained by pairwise 
rankings, a minimum being the lesser of the rank sum 7’, and its conju- 
gate, T’ = (2n + 1)n — T; minimum treatment is the treatment for 
which the rank sum is minimum. Ties were assigned their average 
rank. This gives a multivariate criterion with 15 entries. 


Comparison H, Sf H, Sb H,C L, Sf H,M E¢ 
Minimum Rank Sum 56 57 58 60 62 644 
Minimum Treatment H H H L H L 


Comparison Sb,Sf H,L L,Sb L,M &b,C Sf,M M,C Sb,M Sf,C 
Minimum 

Rank Sum 71 75 75 . 754 80 82 84 100 103 
Minimum 

Treatment Sb H L L Sb Sf M Sb Sf 


From Table 2, a rank sum of 67 is significant at the 5 percent level, 
k = 6, n = 10; hence six comparisons are declared significant. 

The device, used with multiple comparisons procedures, of under- 
lining treatments which cannot be distinguished by their means may 
be adapted to apply to rank sum procedures. Thus, from the test, 
it appears that H and L, and L, Sb, M, Sf and C form two groups 
as a first step; since L can be distinguished from Sf and C, the latter 


group becomes two, namely L, Sb and M, and Sb, M, Sf and C. We 
have: 


H L S$ M Sf C 


Ordering of Sb and M, and of Sf and C was done on the basis of rank 
sums for these paired comparisons though this does not imply that 
this is the only, or even the best, method. 

This procedure is an analogue of Tukey’s w-procedure [19]. The 
significance level is for an experimentwise error rate. Use of Wilcoxon’s 
[21] two-sample test, with its per-comparison error rate, calls for a 
significant rank sum of 78; this will result in four more comparisons 
being declared significant. 


1 


544 BIOMETRICS, DECEMBER 1961 


TABLE 2 
PERCENTAGE PoINTs OF THE MINIMUM RANK SuM 
(AN APPROXIMATION ) 


Number k = number of treatments being tested 
in a 

treatment 3 4 5 6 f 8 9 10 

: 10 7 6 — — — 
10 26 24 23 22 22 21 21 _— 
6 05 “4 23 2 2 — — 
10 36 34 33 32 31 30 30 8 29 
7 05 34 32 31 30 29 28 28 -= 
4 r 10 48 46 44 43 42 41 40 40 
wi 8 05 45 43 42 40 40 39 38 38 
a 10 62 59 57 56 55 54 53 52 
9 05 59 56 54 53 52 51 50 49 
01 52 50 48 47 46 45 — _- 
10 77 74 72 70 69 68 67 66 
10 .05 74 71 68 67 66 64 63 63 
01 66 63 62 60 59 58 57 56 


It is also possible to propose and carry out a sequential procedure, 
an analogue of the Newman-Keuls procedure, which uses several rank 
sums for testing. For this procedure, the above analysis is the first 
step and has separated the treatments into three groups. To proceed, 
we assume that declared differences are indeed real. Hence to test 

“ H versus L, the first group, we may use Wilcoxon’s [21] two-sample test, 
H and L are declared significantly different and the line beneath them 
may be removed. 

Further, compare L, Sb and M using the critical value for k = 3, 
n = 10, namely 74. J versus Sb and L versus M at 75 and 753 are be- 


i 
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yond the 10 percent point but are not quite significant. We cannot 
distinguish among these three treatments. 

Finally, consider the group composed of treatments Sb, M, Sf and C. 
The critical value is 71, k = 4,n = 10. The treatments Sb and Sf can 
be distinguished and we must, then, change the order of Sf and C from 
that proposed when the fixed rank sum test was used. No further 

‘differences will be declared significant by this procedure. We have: 


H L S& M C Sf 


For the Tukey and Newman-Keuls parametric procedures, we find . 
means of 160.2(H), 211.0(L), 267.4(Sb), 278:8(M), 323.2(Sf) and 
326.8(C). Also s, = 18.03 and significant ranges are 51.2, 61.5, 67.6, 
71.9 and 75.2 fork = 2, --+ , 6 respectively. 

For Tukey’s test, the fixed rank sum is 75.2. We find: 


H L M Sf 


For the Newman-Keuls procedure, we find: 
H L M Sf @ 


We now compare the results obtained from applying the parametric 
and non-parametric procedures. 

The fixed rank sum test and Tukey’s test lead to the same con- 
clusions. Both are based on experimentwise error rates. 

Conclusions drawn from the multiple rank sum test and the Newman- 
Keuls test differ as follows. L versus M is significant by the Newman- 
Keuls procedure only; Sb versus Sf is significant by the rank sum test 
only. Otherwise, the procedures lead to the same conclusions. Fifteen 
paired tests have been made. Since L versus M is nearly significant 
and Sb versus Sf is just significant by the multiple rank sum test, it 
would appear that the two methods lead to conlcusions, for this ex- 
ample, that differ only slightly. 

The other multiple rank sum test to be considered is an analogue of 
Duncan's new multiple range test. Table 3 provides critical values. 
A rank sum of 75 is significant at the 5 percent level, k = 6, n = 10; 
hence nine comparisons are significant. Tentatively, we have: 


H L M _ Sf 


Unfortunately, this includes an anomaly since Sb versus Sf is also de- 
clared significant. The same would be true if the Wilcoxon-Mann- 
Whitney test were being used at a significance level (between 5 percent 
and 1 percent) calling for a critical value of 75. 


a j 
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TABLE 3 


PERCENTAGE PoINnTs OF THE MINIMUM RANK SuM FOR DUNCAN ANALOGUE 
(AN APPROXIMATION) 


Number k = number of treatments being tested 
in a 

treatment 3 4 5 6 7 8 9 10 
.10 6100 (1020S 
4 .05 10 10 10 10 10 10 10 10 
10 18 18 17 17 17 17 17 17 
5 .05 17 17 16 16 16 16 16 16 
01 56 6 b—- =— 
.10 27 26 26 26 26 25 25 25 
6 .05 26 25 25 25 25 24 24 24 
.O1 23 22 22 22 22 21 21 21 
.10 37 37 37 36 36 36 36 36 
7 .05 36 35 35 ° 35 34 34 34 34 
.01 32 32 31 31 30 30 30 30 
-10 50 49 49 49 48 48 48 48 
8 .05 48 47 47 46 46 46 46 45 
01 48 42 42 41 #41 #41 «41 «40 
.10 65 64 63 63 62 62 62 62 
9 .05 62 61 60 60 60 59 59 59 
.01 56 55 54 54 53 53 53 52 
.10 81 80 79 79 78 78 78 77 
10 .05 77 76 76 75 75 74 74 74 
.01 70 69 68 68 67 67 67 66 


On the basis that treatments declared significantly different are 
indeed so, we proceed to test treatments M, Sb, C and Sf using k = 4 
and L and M using k = 2. The final result is: 


H 


L 


Sb 


M 


C 


Sf 


This result is the same as that obtained using the Wilcoxon-Mann- 


Whitney test. 


Using Duncan’s (parametric) new multiple range test, we find: 
Sb 


H 


L 


M 


Sf 


Cc 
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Conclusions from Duncan’s test and its rank sum analogue differ 
only in that the parametric test finds Sb versus C to be significant. 

The multiple rank sum test can be adapted to apply to testing 
treatments versus control as well. 


5. APPENDIX—THEORY 


The problem is concerned with pairwise testing of treatments in the 
one-way classification or completely random design. 

The test criteria are rank sums, computed as for Wilcoxon’s [21] 
two sample test, for appropriate pairs of treatments. Rank sums will 
be referred to Tables 2 or 3 for testing rather than to White’s [20] table 
for the Wilcoxon-Mann-Whitney test. 

Two tests will be considered: 

1. The all pairs test, in detail. 

2. Treatments against control, rather briefly. 

For the all pairs test, let X; , 7 = 1, --+ , k be random variables 
measuring some characteristic for each of k samples or treatments. 
Let there be n; observations on the 7-th treatment. Computation of 
the test criterion requires us to: 

1. Rank the X,’s and the X,’s, all ¢ < j, assigning rank 1 to the 

least observation. 

2. Add ranks for the variable with fewer observations to give 7; . 

(There is no loss of generality if we assume n, < m2 < --- < m). 

3. Compute the conjugate of T,; , namely T/;=(n,+n;+1)n,—Ty; . 

The conjugate is the rank total that would be obtained if rank 1 
were assigned to the highest observation. Conjugates are required for 
two-tailed tests. 


Consider (7,2, Tie, Tos , , Te-1,2), Criterion with 


components. Rank tests are based on the assumption that, under H, , 
all permutations of the >> n; observations are equally likely. Hence, we 
must know the number of ways in which (7); , --- , T's-1,2) can be 
obtained. For this, a recursion formula is given in [14]. This provides a 
method, though tedious, of deriving the distribution of (713, --+ , Ts-1,s)- 


The distribution of (T,. , --- , Tx-1,.) has been shown to have the 
following parameters [14]: 


E(T ii) = wis = +0; + 

— wis)” = = nnn, +n; + 1/12, 

— — Bas) = = = on; , 
E(T 


— pa i; — wis) = = —nnn,;/12, 


| 
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ET, — — wii) = = 9, 
Pinay = Pirin = Pains = +7, + + 0; + 1), 
= 0. 
The determinant of the variance-covariance matrix is: 


(II ni (>> n; + yy12, 


The elements in the inverse of the variance-covariance matrix are: 
Corresponding to the variance of 7; , 
12(>> nme n;)/nn,( >> N. + 1). 
Corresponding to a covariance with the 7’s having a common 
subscript, h, in the same position: 
—12/n,( + 1). 
Corresponding to a covariance with the 7’s having a common sub- 
script, 7, in different positions, for example 7; and T7;; : 
12/n( + 1). 
Finally, the element corresponding to ¢,,.;; is zero. 
The determinant of the variance-covariance matrix may be evaluated 
as follows: from the ith row of the determinant, factor n,n,.;/12, 
t= 1,--- ,k — 1; from the ([k — 1] + 2)-th row, factor n.n-,,;/12, 
1,---,k — 2; ; from the last row, factor n,_,n,/12. The 


k 
product of these factors is [J n‘~'/ 12°), 
The entries in the determinant which is the other factor may be 
described somewhat crudely as follows. 


The i-th diagonal block, i = 1, --- , k — 1, which contains the 
variances and covariances of the T7';;’s, fixed i and j > 7, will be 


Nn Ni+2 Ny 


Nest tl Ne 


This is a (k — i) X (k — i) block. 
The block consisting of the same rows and the first (k — 1) columns is 


F — 2 columns k — 7 columns 
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The next block to the right consists of k — 2 columns and is: 


--- m 00 --- 0 | 

0 --- 0 
— 3 columns k—decolumns | 


The pattern is now clear. 
The first block to the right of the i-th diagonal block consists of the _ 
next (k — [¢ + 1]) columns and is 


| —Ny43 eee 

te 

0 
_ 

The block to the right of this is 

“Ries “Mies *** 

0 Nise O 


Again, the pattern is clear. 

At the next step in evaluating the determinant, from column 1 sub- 
tract columns ({k—1]+1) through ((k—1]+[k—2]). To column 2, add 
column ({k — 1] + 1) and subtract columns ((k — 1] + [k — 2] + 1) 
through ((k — 1] + [k — 2] + [k — 3]). In general, to column 
i,i = 1, --- ,k — 1, add the z — 1 columns described as those con- 
taining the first column of the diagonal blocks previously set out, and 
subtract the sum of the next (k — [7 + 1]) columns, that is, the columns 
including the (7 + 1)-st diagonal block. 

The result of the above operations is that the first diagonal block 
has >> n; + 1 as common diagonal element and zeros elsewhere. The 
other blocks containing the first k — 1 columns are of the form previously 
described but with —n, and n, replaced by —(> 7, + 1) and p> nti 
respectively. 

At the next and final step, leave rows 1, --- , & — 1 unaltered. 


: 
per 
oy, 
te 
4 
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Consider the next k — 2 rows which we now call the first set, the next 
k — 3 rows called the second set, and so on, the i-th set consisting of 
(k — [i + 1)) rows, beginning with row >-‘., (k — a) + 1. To the 
j-the row of the i-th set, add row 7 and subtract row 7 + j, these rows 
coming from the first k — 1 rows of the determinant. 

The result of these operations is a determinant with >> n, + 1 in 
the first k — 1 principal diagonal positions, ones elsewhere in the princi- 
pal diagonal, and zeros below the principal diagonal. Hence the de- 
terminant is as given. 

That the inverse elements are correctly given may be checked by 
multiplying the matrix by the given inverse. 

For the treatments versus control procedure, let X;,7 = 0,1, --- ,k 
be random variables measuring some characteristic of a control and 
k treatments with n, observations in the i-th sample. 

Computation of the test criterion requires us to: 

1. Rank jointly the X,’s and X,’s, i fixed, giving rank 1 to the 

least observation. 

2. Add ranks for the variable with fewer observations, here assumed 

to be the check, to give T; . 

3. Compute the conjugate, 7% . 

Consider (7, , --- , 7,). Again, a recursion formula for finding the 
number of permutations which give rise to a specific value of 
(T; , -+- , Ty) is given in [13]. 

The distribution of (7, , --- , T',) has been shown to have the follow- 
ing parameters [13]: 


E(T,) = = + 2; + 1)/2, 
KT, == NN (NM +n; + 1)/12, 
E(T; — — wi) = = nonn,;/12, 
pis = + + +0; + 1). 
It may also be shown that the determinant of the variance-co- 


variance matrix is 
+ ny + 1)12. 
The diagonal and off-diagonal elements of the inverse are, respectively, 


12( + 1) / [ mane + Na + 


[nine + n; + 1) |. 


and 


= 
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To evaluate the determinant of the variance-covariance matrix, 
factor non;/12 from the i-th row of the determinant. 
This gives: 


k 
k-1 N +n +1 Ne 
No IIx. 
0 


ny Ns No tm + 1- 
Next, subtract the k-th row from the i-th row, i = 1, --- ,k — 1. 
We obtain 


nt +1 0 + 1) 
Ne Ntmtil 
Finally, obtain a new k-th column as the sum of all columns. It is 
then apparent that the given determinant is correct. 
That the given inverse elements are correct is seen by multiplying 
the matrix by the stated inverse. 


From the above information, it is possible to tabulate probabilities 
for the parent distribution and derived distributions of interest. 
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THE ESTIMATION OF REPEATABILITY AND 
HERITABILITY FROM 
RECORDS SUBJECT TO CULLING 


R. N. Curnow! 


Agricultural Research Council Unit of Statistics, 
University of Aberdeen, Aberdeen, Scotland? 


1. INTRODUCTION 


In any animal breeding selection programme, estimates of repeat- 
ability and heritability are needed to choose between the various 
selection schemes available and also to ensure that the highest possible 
genetic gains are obtained from the chosen scheme. Estimates of 
repeatability and heritability are generally subject to large sampling 
errors. Therefore, the most efficient methods of estimation should 
be used even if they do involve rather lengthy computations. In this 
oaper, the maximum likelihood estimation of repeatability and herit- 
ability from records subject to culling will be considered. The more 
usual regression estimators are often very inefficient compared with 
these maximum likelihood estimators. 

Suppose that we wish to estimate the repeatability of lactation yield 
in a herd of dairy cattle. We shall assume that only first and second 
lactation yields are available and that, if there had been no culling 
(t.e., if all the cows had had second records as well as first records), 
the first and second records would have been normally distributed over 
the herd with means py, and yp, , variances o; and a; and covariance 
between the two records of the same cow o,,. The first and second 
records of cow 7 will be written y,;, and y;, respectively. The assumption 
of normality for the distributions will probably be a reasonably good 
approximation unless the herd can be split into groups so that any two 
cows in the same group are much more alike than two cows in different 
groups. These groups may, for example, be groups of daughters of the 
same sire or groups of cows according to the year in which they gave 
their first record or the month in which they calved. Methods are 
available for making allowances for such groupings, but they will be 


assumed absent in the rest of this paper. Very rarely will the culling 
1The work on this paper was undertaken while the author was visiting the United States on a 
Harkness Fellowship of the Commonwealth Fund. 
2Pregent address: Unit of Biometry, University of Reading. 
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in the herd be sufficiently intense or sufficiently highly correlated with 
future milk yielding capacity to affect seriously the normality of the 
distributions. 

Repeatability is defined as the correlation between two different 
records of the same cow and is therefore 


P= 12/0102 


Since we are considering only first and second records, the question 
of whether the repeatability is the same for all pairs of records will not 
be discussed. The assumption is frequently made that oj = o} and, 
therefore, that repeatability is the same quantity as the regression 
coefficient of second records on first, 82; = o12/0; . The assumption 
is made, for example, in the formula given by Lush [1945] for comparing 
cows with a differing number of records and in the formula given by 
Lerner [1958] for the ratio of the heritability of the mean of n records 
to the heritability of a single record. It was also used by Henderson, 
Kempthorne, Searle and von Krosigk [1959] in their discussion of the 
disentanglement of environmental and genetic trends from records 
subject to culling. When of ~ o} , an estimate of 8,, is needed for 
prediction purposes but for other purposes estimates of oj and o2 may 
also be required. We shall assume in this paper that 07 = 03 = o’, 
say, and therefore that 6., = p. A logarithmic transformation applied 
to the data may sometimes be useful in satisfying this assumption. 

p is generally estimated by b, the sample regression coefficient of 
second records on first. There are two reasons for this. First, b is 
very simple to calculate and, second, it is an unbiassed estimator of p 
despite any culling that may have been based on the first records. 
However, when o; = o} and the first records of all cows are available, 
whether or not they have second records, b is not the maximum likeli- 
hood estimator of p. The variance of the second records about their 
regression on the first estimates o*(1 — p’) and the variance of all the 
first records estimates o*. These two estimators can be combined to 
give an estimator of p’. Maximum likelihood makes use of this in- 
formation as well as the information given by b. 

We shall derive the maximum likelihood estimators of p and a’. 
We shall show that the efficiency of b as an estimator of p, relative to 
the maximum likelihood estimator, is fairly low for values of the various 
parameters that may well occur in practice. The bias of the maximum 
likelihood estimator is shown to be small. This suggests that the 
maximum likelihood estimator may be worth calculating. The com- 
putations involve only the solution of a cubic equation. A section is 
devoted to an approximate check of the assumption that o? = o? . 
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Attention has been confined so far to the estimation of repeat- 
ability. The methods to be discussed could also be applied to the 
estiination of heritability from parent-offspring records. The assump- 
tion oj = o; means that, had there been no selection of parents, the 
variance of the parents and of the offspring would have been equal. 
In heritability studies, the maximum likelihood method makes use of 
information about animals that are not parents of animals which also 
have records. In milk yield studies, this would include information on 
dams that have only male calves. 


2. THE MAXIMUM LIKELIHOOD ESTIMATION OF REPEATABILITY 


All the cows have first lactation records but they do not all have 
second lactation records. We shall assume that the probability that a 
cow has a second record depends on its first record but not on any other 
character correlated with the second record. This rules out culling 
based on information about relatives or on characters such as per- 
centage butter-fat. However, the effect of such culling on the estimates 
will often be very small and could probably be safely ignored. 

t N cows have a first record and n < N cows have a second record. 
Numbering the cows with a second record i = 1, 2, --- , n and the 
cows without a second record i = n + 1,n + 2, --- , N, the first records 
can be written 


Yar = 1,2,-:- ,N) 


and the second records . 
Yi2 = 1,2,--+ ,n). 
The N first records are normally distributed with mean y, and variance 


a. Because the culling is based only on first records, the distribution 
of a second record y,2 , given the first y,;, , is independent of the dis- 
tribution of the first record and is normal with mean yp, + p(yir — #1) 


and variance «(1 — p’). The likelihood of all the records is therefore 


~ 20° 
1 [yi2 — we — — 


and the log likelihood is, apart from a constant, 


i 
4 
4 
N 
L = 
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InL = no —p) 


20° a — p) (2.1) 


Kempthorne and von Krosigk (Henderson, Kempthorne, Searle and 
von Krosigk [1959]) suggested the use of maximum likelihood to 
estimate repeatability from records subject to culling. They derived 
the log likelihood when the cows were classified into groups and when 
the number of records for each cow was perfectly general. They derived 
very lengthy equations from which the maximum likelihood estimates 
of the parameters could be determined. In this paper, we are studying 
in greater detail the special case when there is no grouping of the cows 
and when there are at most only two records per cow. 
(2.1) can be written 


n(p + 23 — 2 
2o"(1 — (2.2) 
20° 20°(1 — p’) 


Ya — 9)” > Ya — 


The maximum likelihood estimators of 4, and yw, are clearly given by 


=9 and = — — 


Also 
ainL)_ _@+N) + Ns? n(p + — 22,2) 
2¢ 20°(1 — p’) 


(2.3) 


N@ — — we — eG — 
20°(1 — 


4 
i { 
where 
= 
22 ay n “12 n 
N n 
Ya Yar Yi2 
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and 
Ain L)_ __ mp + 22) — (1 + 
dp o(1 — p) 


nplGo — — — 
o (1 p) 
By setting (2.3) and (2.4) equal to zero and substituting the maxi- 
mum likelihood estimators of u, and yp, , A satisfies the cubic 


(Ns* — — (N 
+ [(n + — Ns’ — — (n+ = O. 


This equation will have one and only one root between --1 and +1 
of the same sign as the simple estimator b = 2,,/2;. There may 


be two other roots between —1 and +1 of opposite sign to b. The 
formula for é” is 


(2.5) 


A(l — 


When there is no culling, s? = 2? and 


2219 


When of = o; , the maximum likelihood estimator of the correlation 
coefficient has the arithmetic rather than the geometric mean of the 
two sample variances in the denominator. 

To obtain the asymptotic variance of 6, we need the expected values 
of the second-order partial derivatives of In Z with respect to wu: , us, p 
and o”. These can be derived from a knowledge of the expected values 
of s’, 2} and 2... 2? depends on the method of culling and so will 
be taken as fixed. For reasons to be given later, we shall calculate the 
expected values when o? ¥ Clearly, E(s’) = (VN — 1/N)oi. By 
writing 


» 


Yi2 = Me + to V1I- p 
1 
so that e;., has a standard normal distribution independent of yi. , 
n 2 
1 


and, therefore, for given 27 , 
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Similarly, 


Also, with of = 0; , E(g — = o°/N, — us — — = 
o°(1 — p’)/n and E{(g, — — — — m))} = 0. After con- 
siderable algebra the asymptotic variance—covariance matrix for A, , 
fA. , # and @ is found to be V, where 


n n 0 nH. — wn) 


~o(l — o(1 — p) o(1 — p’) 
n+wN np 


n 2 2 ny, i)” 
— + + me 


The asymptotic variance of f is 


= (N_+n)(1 — 

n{(N + n)(1 — + 2Np'] 
Writing n/N = S, so that (1 — S) is the culling intensity, and 
x?/c” = 1/c, where c measures the effect of the culling on the variance 


of the first records and is therefore a possible measure of the efficiency 
of the culling, 


‘| 
We shall now derive the approximate bias in # as an estimator of p. 
By substituting 6 = p + 64 in (2.5) and ignoring terms in 3 and 3’, 


3(Ns’ —nz,”)p —2(N —n) 2,2 + [(n + 23 )] 


(2.7) 


(2.8) 
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Approximating #(6) by taking expected values separately in the numer- 
ator and denominator, 


+ — #) + — N) + + N — 3n)] 


To order 1/n, 


E(s) = = 
+ 29 | 
Therefore, from (2.7), the approximate ratio of the bias to the standard 
deviation is 
E(é)_ p(l — S) 


1 | 
If n is reasonably large, the bias in / is unlikely to be serious. The bias 


will become relatively more important when estimates of p from different 
sources are pooled. 


(2.9) 


3. ESTIMATION OF o? AND WHEN o? 03 


In the previous section we derived the expected values of 8’, 23 
and for given 27 when o? They were 


E(s’) = 


E(2,2) = = 


Apart from the usual slight differences in the multipliers, important 
only for small n and N, the following estimators of of , 8 = po2/c, and 
o3(1 — are the maximum likelihood estimators when , 


= Ns’*/(N — 1), 
Bar =b= 212/24 


{ 
‘ 
4 
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n 
— p) = [23 — 22/2). 


From the expected values of .s’, 23 and 2,2 given in the previous 
section these estimators are all unbiassed. The maximum likelihood 
estimator of is 


i 2 2 2 2 
(3.1) 
No exact method is available for constructing a confidence interval for 
o;/o; except when n = N, i.e., when there is no culling (see Curnow 


[1957] for references and for a method sad Sa when the data are 


grouped). The asymptotic variance of (03/c) could be derived and 
used to provide an approximate confidence interval. This would give 
some indication of the importance of the assumption that of = o} . 
A simpler, but much less sensitive, test of whether o> > o{ could be 
based on the fact that the quantity 


N(n — — p) 
which compares a x’-value from the variation of the second records 
about their regression on the first with a x’-value from the variation 


of all the first records, has an F-distribution with n — 2 and N — 1 
degrees of freedom. 


(N — 1)n 22 — 22/21 
N(n — 2) 


significantly greater than F = 1 would suggest that o3(1 — p’) > a3 
and, therefore, that > . 


4. THE EFFICIENCY OF THE SIMPLE REGRESSION ESTIMATOR 
OF REPEATABILITY 


The statistic b = 2,./Z{ is very easy to calculate and is the one 
generally used to estimate repeatability. It is the only possible esti- 
mator if first records are available only for those cows also having second 
records. Providing that culling is based solely on the first lactation 
yields, b is always an unbiassed estimator of B., = po2/o, , but an un- 
biassed estimator of p only when oj = o;. In this paper we are assuming 
a; = o,. The variance of b is 


V(b) = o°/ndi = c/n. (4.1) 
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From (2.7), the asymptotic efficiency of b relative to the maximum 
likelihood estimator, #, is therefore 


The values of this efficiency are shown in Table 1 for various values 
of S, p and c. 

The following considerations suggest that c is unlikely to be greater 
than c = 2. Let the N cows be reduced to LN by accidental factors 
uncorrelated with the level of yield. Then a proportion P=n/LN=S/L 
can be selected on the basis of first lactation yield. Assume that the 
selection is of the proportion S/L of the herd having the highest yields 
and that n and LN are sufficiently large that the effect of the selection 


TABLE 1 


Tue Asymptotic EFFICIENCY OF THE SIMPLE REGRESSION EsTIMATOR OF 
REPEATABILITY RELATIVE TO THE MAXIMUM LIKELIHOOD EsTIMATOR 


Overall Selec- c = 09/3? 
tion Intensity Repeatability 
(S = n/N) (p) 1/2 


4/3 


1/4 0 1.00 1.00 1.00 . 1.00 
1/4 0.95 0.90 0.88 0.82 

1/2 0.79 0.65 0.58 0.48 

3/4 0.49 0.33 0.27 0.20 

0 0 0 0 


1/2 0 1.00 1.00 1.00 1.00 1.00 
1/4 0.96 0.92 0.89 0.85 0.74 

1/2 0.82 069 O68 0.53 0.36 

3/4 0.54 0.387 40.30 0.23 0.13 

0 0 0 0 0 


3/4 0 1.00 1.00 1.00 1.00 1.00 
1/4 0.96 0.93 0.91 0.87 0.77 
1/2 0.84 0.72 0.66 0.57 0.40 
3/4 0.58 040 0.384 +O.25 0.15 
0 0 0 0 0 


1 0 1.00 1.00 1.00 1.00 1.00 
1/4 0.97 0.94 0.92 0.88 0.79 
1/2 0.86 0.75 0.69 0.60 0.43 
3/4 0.61 0.44 0.37 0.28 0.16 
0 0 0 0 0 


“Be 
i 
| 
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can be approximated by the truncation of the top P proportion of an 
infinite normal population with mean yz, and variance o”*. The variance 
of the selected first records will then be 


z= — rv — T)], 


where vy = Z/P and Z and 7 are the ordinate and abscissa at the point 
above which lies a proportion P of the population (I’inney [1957] derives 
this formula and provides a table of values of »’ = »(v — T) for various 
values of P). The division of the overall selection intensity, S, between 
L and P is never complete and is certainly never known exactly. How- 
ever, P is generally considerably larger than S, 7.e., there is a large 
amount of random culling and not so much actual selection. The 
value of c is 
: 2 


For various P, it takes the following values: 
P=0.95 09 08 0.1, 
c= 1.24 1.40 1.72 2.03 2.387 2.75 5.91. 


In practice, P is unlikely to be less than P = 0.7 and therefore c is 
unlikely to exceed c = 2. Table 1 does include c = 4 to illustrate the 
effect of intense selection based on yield and c = } to illustrate the 
effect of a selection scheme that results in an increased rather than a 
decreased variance. The efficiency is the same for p as for —p and so p 
is shown as positive in the table. In practice, p is very unlikely to be 
negative. 

When p = 0, the efficiency of the simple regression estimator is 1. 
When p = 1, it is 0. For intermediate values of p, the efficiency in- 
creases rather slowly with S for fixed c, but decreases fairly rapidly 
with increasing c. For all values of c and S, the efficiency is very low for 
high values of p. As an example of the efficiency likely to occur in the 
estimation of the repeatability of lactation records, the efficiency is 
0.66 when c = $, p = 3. and S = 3. In this case, the maximum likeli- 
hood estimator would certainly be worth calculating. 


5. THE ESTIMATION OF REPEATABILITY BY AN ANALYSIS 
OF VARIANCE 


Sometimes the effect of culling based on the first records is ignored 
and p estimated from a least squares analysis of variance of the N + n 
first and second records. Table 2 shows this analysis of variance. 
The “cows and periods” sum of squares has been split into both “cows 
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adjusting for periods and periods ignoring cows’ and “cows ignoring 
periods and periods adjusting for cows’. By a period difference is 
meant an average difference between fitst and second lactation records. 
The expected values of four important mean squares are shown. A, , 
A, and A; are defined at the foot of the table and vanish only if the 
culling is random with respect to the first records. Two cases need to be 
distinguished. In the first we assume u, = mu, and in the second we do 
not. We shall use the reference numbers at the right of the table. If 
#; = M2 , Mean square 2 and a mean square obtained by pooling sums 
of squares 3 and 4 are used to estimate o” and p. This is exactly equiva- 
lent to an analysis of a one way classification with unequal sub-class 
numbers (Snedecor, [1956], §10.16). If u, + uw, , mean squares 1 and 4 
are used to estimate o” and p. There seems to be little justification 
for using mean square 2 with mean square 4. One or other is biassed 
or inefficient according as OF = by. 

The bias in any of the above methods is difficult to determine without 
some knowledge of the effect of culling on the expected values of 
— [(n — 1)/n]o’, — — — o*/nand — o*/N. 
However, these unknown biasses are clearly undesirable. Since the 
method of estimation described in this paper is the maximum likelihood 
method whether or not the culling is random with respect to first 
records, it is to be preferred. If culling is random with respect to first 
records, c can be substituted as c = 1 in the formulae for the asymptotic 
bias and variance of the maximum likelihood estimator and the asymp- 
totic efficiency of the regression method. 

Wadell [1961] has given a method for estimating p when the culling 
is equivalent to truncation of the first lactation yields. His estimate 
of p is a function of g, , J. , s°, 2? and the mean squares 2, 3, and 4 of 
Table 2. The estimate has a negligible bias for large values of n but 
its variance is not given and therefore its efficiency is unknown. 


6. A NUMERICAL EXAMPLE 


Dr. A. E. Freeman of the Department of Animal Husbandry, Iowa 
State University has kindly made available to me some Complete Herd 
Improvement Registry Records of an Iowa Board of Control herd of 
Holstein-Friesian cattle. The records included some first and second. 
lactation yields made in various years but expressed as deviations trom 
the herd average for each particular year. The cows were milked 
twice daily and their yields expressed as mature equivalent 305-day 
lactation yields. 

The values of the various relevant statistics for a sample of the 
records were: 


BS 
‘ 
4 


BIOMETRICS, DECEMBER 1961 


TABLE 2 
ANALYSIS oF VARIANCE OF First anp Seconp Recorps 


Expected Values of Mean Squares 


Sums of Squares 


(22 + 222 — + Ne 


(n —(N + 3n) 
22. — 2) + NE 

N 
V+ n @ — — 92) = 
+ 


5 — p) +: (u2 — +n(1 — p) As 


+2 1 — 


Adj. for (N. (Ng + 
mean N+n 


4, = 7, 


1 


N 2N 
+ + — D (Dg — 0° /N), 


4, = 1 p)(D*® = o°/n) — D and D= 


N=220, n=150, S = 0.682; 
g=-1.118, = 1.147, = —0.207; 
* = 286.55, = 270.51, 125.85 and = 354.98. 
The maximum likelihood estimator of o3/a7 (3.1) is (02/c?) = 1.25. 
This value is almost certainly not significantly greater than o2/o? = 
However, it is sufficiently large to raise some Coubts about the assump- 
tion that ¢, = o| . The estimate of o3(1 — p’) is greater, but not 


8 


— 
| 
nt+N-2), 
Periods 
Cows and N 5 + 22,. + Ns’ | 
— p) 
Cows ign. 
| Periods | 2 
|! 
Periods 
adj. Cows | : 
4 | 
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significantly greater, than the estimate of «7. However, for illustrative 
purposes we shall assume o3 = o} . 
Equation (2.5) for / is 
p° — 0.39224” + 5.12576 — 2.0728 = 0. 
The simple regression estimator of p is 
b = 22/2; = 0.465. 


Three cycles of an iteration based on formulae (2.8) and with 6. = b 
as the initial solution show that the value of 6 is, to three decimal 
places, 

6 = 0.404. 


This is the only real root of the maximum likelihood equation for #. 
The two estimators together with their estimated standard errors 
[(2.7) and (4.1)] are therefore: 


b = 0.465 + 0.084 


and 


6 = 0.404 + 0.069. 


The latter standard error is an asymptotic standard error. The esti- 
mated efficiency of the regression estimator is only 67 percent and so 
é is probably well worth calculating. The estimated asymptotic bias 
(2.9) is very small, 

E(6) = 0.0004. 


The maximum likelihood estimator of o” (2.6) is ¢ = 314.46. This 


agrees well with the value s’ = 286.55 used above in estimating E(6) 
and the standard errors of b and #. 
The estimated value of c is 


é = 3°/Zi = 1.059. 
This value of ¢ suggests that there has been very little selection based 


on first lactation yield in this particular sample of records. The 
standardized selection differential 


— = 0.134 
suggests that the effective selection intensity P was near 0.93. 


7. SUMMARY 


The maximum likelihood estimation of repeatability from first and 
second lactation records of a herd subject to culling is discussed. The 
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method is applicable only when it can be assumed that, had there been 
no culling, the variances of the first and the second lactation records 
would have been equal, and when all the first records are available 
whether or not the cow has a second record. The efficiency of the more 
usual estimate of repeatability, based on the regression of second records 
on first records, is shown to be low when the repeatability is high. In 
many cases the calculation of the maximum likelihood estimate would 
seem to be well worthwhile. The use of the method in estimating 
heritability is mentioned. An illustrative example is given. 
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THREE CLASSES OF UNIVARIATE DISCRETE 
DISTRIBUTIONS 


C. G. Kuarri I. R. Pate. 
M. S. University of Baroda, Baroda, India. 


1. INTRODUCTION 


Families of descrete distributions have been developed and studied 
by many authors, including, Neyman [1939], Feller [1943], Skellam 
[1952], Beall and Rescia [1953) and Gurland [1957, 1958]. These 
families are of three types: 


Type A: ga(z) = exp {h@)}, 
Type B: ga(z) = {h@)}", 
Type C: gce(z) = c log {h@)}, 


where g(z) represents a probability generating function (p.g.f.) and 
h(z) is a p.g.f., except possibly for additive and multiplicative constants. 
The aim of this paper is to set up formulae for certain statistics for 
these types. It is hoped that these will be of use to reseach workers in 
practical fields, who will be formulating compound and generalised 
distributions of these types by using specific forms of h(z). 


2. RECURRENCE RELATIONS FOR PROBABILITIES 


Notations: Let the r-th derivation of f(z) be denoted as f‘”(z). 
Also let 


and 


a, = {h'(2)/r!} leno - (2) 


By the definition of g(z), it is clear that P, denotes the probability 
of the r-th count in g(z) and z, , the probability of the r-th count, 
excepting possibly for additive and multiplicative constants. 

Type A: Here 


ga(z) = exp {h(z)}. 


1 Present address: Statistical Officer (I), Rajkot, India. 
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= 


Successive differentiation leads to 
k=1 
Hence, on letting z = 0, we have 
P, = with Po = exp (m). 
k=1 


Type B: Here 
gu(2) = {h(z)}”. 


Successive differentiation of h(z)g$?(z) = nga(z)h"’(z) leads to 


k=1 


Hence, on letting z = 0, we have 
P, = (nk with Py = 7. 
k=1 


Type C: Here 


gc(2) = c log {h(z)} where c = log {h(1)}”’. 
Successive differentiation of h(z)g‘(z) = ch (z) leads to 


k=1 


Hence; on letting z = 0, we have 


r-1 
P, = {rer, — — with = clog m. 


3. FACTORIAL CUMULANTS 
Notations: Let 


(@} 


Mir 
Mi, 


and 


= {(d/dz)" log g)} |.-: - 


I'rom the definition of factorial cumulants, it is clear that yu/,, is the 
r-th factorial moment of h(z) if h(z) is a p.g.f. and M{,, and K,,, are 
respectively r-th factorial moment and r-th factorial cumulant of g(z). 
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(3) 


(4) 


(6) 


(7) 


(8) 


(9) 
(10) 


(11) 


k=1 
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Type A: Here log {g4(z)} = h(z) and hence, 
Ki, = Bir) 
Type B: Here 


o(z) = log ga(z) = n log h(z). 


Hence, on using (11), it is clear that r-th factorial cumulant of g,(z) 
is n times the r-th factorial cumulant of h(z) if h(z) is a p.g.f. Also on 
using (7) (with necessary modifications), (9) and (11), we have 


r-1 


~ ] 

Type C: Here, it is easy to give a recurrence relation for factorial 
moments M/{,, rather than factorial cumulants. By using (7), (9) and 
(10), this relation can be shown to be 


(14) 


where uf, = A(1). The factorial cumulants can be obtained from (14) 
by using the relations 


Kis; = — + , etc. 


4. SPECIAL CASES 


Notation: If a variate has either a Binomial or a Negative Binomial 
law as a special case, we say that it has a general binomial law. 

Type A: From the form of the Poisson-Binomial, Negative-Bi- 
nomial, generalised Polya Aeppli, Beall and Rescia [1953] and Neyman’s 
Types A, B, and C, it is clear that they belong to this form. It is to be 
noted, that, if in the classical problem of egg masses and larvae, the 
egg masses have a Poisson distribution with p.g.f. exp [A(z — 1)] and the 
larvae within an egg mass, a distribution with p.g.f. w(z), then the 
distribution of the larvae over the whole field is exp [A{w(z) — 1}] 
which is of Type A with h(z) = A{w(z) — 1}. 

Some important distributions with their recurrence relations are: 

(i) Poisson-Hypergeometric distribution [1958]: Here 


w(z) = (acme — Bw (15) 


r= 


where 


a.) = aa + (*) =k /rl, 
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k(k — 1) ---(k -—r+1), = 1, = 1 and a,B, m 


and k are such that w‘”(0)/r! = z, is positive. The above distribution 
was first given by Gurland [1958]. The recurrence relation in prob- 
abilities is 
P, sx,P,_./r (16) 
s=1 
where P, = exp {A(m — 1)}, and 
_atB+r —2 — ma —k + 2r — 3) 
r(l — m) 
_ 


Tr-1 


(17) 


For the particular cases of the above distribution, we may refer to 
Gurland [1958]. The recurrence relation for factorial cumulants is 


Ki, = ayn + B) (18) 
(ii) Poisson-Power series distribution: 
wz) = Daz’, (19) 


where a,’s are constants such that w(z) is convergent for some z. The 
recurrence relation for probabilities is 


P, = with P, = exp — 1)}. (20) 


e=1 


The above distribution was first given by Maritz [1952]. 

Type B: Let the distribution of egg masses be a general binomial 
with p.g.f. (1 — p + pz)”. Then the distribution of the larvae over 
the whole field is [1 — p + pw(z)]", ie. the Type B distribution with 
h(z) = {1 — p + pw(z)}. When n = 1, this leads to the distribution 
with p.g.f. w(z) with an addition (or subtraction) of zeros. If n > 1 
or n < 0, this can be regarded as n-th confluent of w(z). 

Some important distributions with their recurrence relations are: 

(i) G. Binomial-Hypergeometric distribution: Here w(z) is the same 
as given in (15). This was first stated by Gurland [1958] in the special 
case when n < 0, p < 0. The recurrence relation in probabilities is 


P,=p (ns — r + 8)2,P,_,/rdo , (21) 


i 
“4 
dr 
dm 
‘a 
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where Py = a} , G@ = 1 — p + pmo and z,’s are defined in (17). 

(ii) G. Binomial-G. Binomial distribution: Here w(z) =(1—m-+mz)*. 
Hence, 7; = m(k — i + 1)x;-,/i(1 — m), with m = (1 — m)*, 
a = 1— p+ pm and Py = aj. The recurrence formula has the same 
form as in (21). 

(iii) G. Binomial-Poisson distribution: Here w(z) = exp {A(z — 1)}. 
Hence, 7; = A‘mo , = exp (—A), @ = 1 p+ pr and P, = 
The recurrence formula has the same form as in (21). Whenn = 1 and 
p{l — exp (—A)} = 86, this was called by A. C. Cohen [1960] an ex- 
tension of a truncated Poisson distribution. 

Type C: Let the distribution of egg masses be a logarithmic dis- 
tribution with p.g.f. log (1 — Az)/log (1 — A). Then the distribution of a 
larva over the whole field will be log {1 — Aww(z)}/log (1 — A), i.e. the 
Type C distribution with h(z) = 1 — \w(z). The important distribu- 
tions are obtained by considering w(z) as hypergeometric function 
G. Binomial, Power-series etc. 


Choice of a distribution under the condition of no migration: 


Let us suppose that the different sites of a colony are distinct and 
countable, and let there be no migration between sites of a colony. 
Assuming the same probabilities of arriving at a particular site by any 


organism, the distribution of r (when r is fixed) organisms in a particular 
site is Binomial. Now let us suppose that one or more organisms 
arriving at the colony follow the truncated Negative Binomial law. 
(This may be true under the wide applicability of the Negs.tive Binomial 
in biological data, e.g., Bliss [1953], Evans [1953]). Then, it is easy to 
show that the p.g.f. of the organisms in a particular site (when there is 
no migration) is 


1 — 6+ Of{(1 + m, — mez) — (1 + m)“}/{1 — (1 + 


where 0 < 6 < 1, m, > 0,k, > O;ie. 1 — p + p(1 + m, — m,z)™ 
for0 < p < {1 — (1+ m,)™}7’, m, > 0 and k, > 0. Now suppose 
that the independent results of n sites are combined together. Then 
the p.g.f. of the distribution of organisms is 


—pt+p(l +m — 
which is a particular case of G. Binomial-G. Binomial. The above 
distribution can be named as Binomial-Negative Binomial. 
5. METHODS OF ESTIMATION FOR G. BINOMIAL-G. BINOMIAL 


Here the p.g.f. is [1 — p + p(1 — m + mz)'J’. 
Method of moments: Let T = Kya R= Ky) Ki1)/Kia and 


. 
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S = Kj,,/Ki2,. The approximate value of n is obtained from 
n(TS — 2k? + R) + nS(TS — 6R + 6) + SR — 2) = 0. 

Then k, m and p are estimated from 

k = {S?+ nS + (2 — R)}/n{n(1 — R) — S}, 


m = (nKj2, + Ki) /(k — 1)nKi, 


and p = K,,,;/knm. (22) 


Method of maximum likelihood when n and k are fixed: 


Similar to Sprott’s results [1958], we have the maximum-likelihood 
equations as 


7 =nkpm and L(m) = (23) 


with = (d/dh)L(m) = a,G(r)[m {1 — kl — 
{1 + (1 — %)[mk(1 — p)]"'h'FA G(r)}], where 7 is the sample mean, 
a, is the frequency at the r-th count, N is the total frequency, 
G(r) = (r + UP,.:/P,7, A G(r) = G(r + 1) — G(r) and m, # are 
maximum likelihood estimates. 

From the first approximate estimates 7p’, m’, the new corrected 


Alt 


values 7”, 7’’ are estimated from 


= mh! — and = F/nkmh"’. (24) 


Sample zero frequency when n and k are fixed: 
Here the estimates p and m are obtained from 


a, = N(l —p+p(l — m)*)” and = nkpm. (25) 


It may be noted that when n = 1, the two equations in (25) are the 
same as those in (23). 


6. EXAMPLE 


In order to illustrate how the above discussion can.help an experi- 
menter in the field of curve fitting, we fit here Binomial-Negative Bi- 
nomial to the data in Distribution 1 of MacGuire et al. [1957, Appendix]. 
From the data, the first four factorial cumulants are K,,; = 2.5900156, 
0.6877630, 0.0218497, Kia == 0.9080044. 

Hence, on taking n = 1, the solution of k, by the method of moments, 
correct to first decimal place is —12.0. Then the various estimates 
of m and p are m = —0.2204963, p = 0.9788582 by the method of zero- 
cell frequency or maximum-likelihood and m = —0.2196584, p = 
0.982591 by the method of moments. 
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The fits are shown in Table 1. The fits of other distributions are 
given along side for reference purposes only. 


TABLE 1. 


Binomial-Neg. Binomial 
Obs. Negative 
fre- by M.L. or | by method | Binomial 

quency zero fr. of moments {1957} 


355 .000 
622 .478 
730.994 
616 .306 
417.545 
241.296 
123 .567 
57.406 
24.632 
9.395 
7.301 


© 


19.184 


2Expected frequencies are from 10 and above 
3Probabilities are calculated from the p.g.f. exp { —(a +b) +az+ bz?} where a = 1.9022526 and 
6 = 0.3438815. 


7. ASYMPTOTIC EFFICIENCIES 


Here we give the asymptotic efficiencies for the various methods 
of estimation for m and p only in a G. Binomial-G. Binomial distribution. 

The determinant of the information matrix up to order N~' for the 
maximum likelihood estimates m and p is 


Dy.» = N’n*[npk{k(1 — p) + (1 — m)m"jR 


— {kl — p) — 


where R = —1 + >> G(r)P, and G(r) = (r + 1)P,4:/nkpmP, . 
The detorminant of the covariance matrix up to order N~* for the 
moment estimates of m and p is 


D = (K,K, + 2Ki — K3)/[Nk(k — 1)nmK,]’, (27) 


(26) 


where 


K, = kn mp, K./K, = + (k 1)m (K,/n), 


Count Poisson; | Poisson 
per Binomial | Power 
plot (1° 37] Series? 
355 346.445 | 324.30 | 341.84 | 339.072 = 
600 628.153 | 660.37 | 644.37 | 645.036 a 
781 | 735.340 | 734.06 | 728.03 | 730.150 “ 
567 618.023 | 610.45 | 609.14 | 610.886 a 
441 417.393 | 408.82 | 415.60 | 416.079 a 
245 240.550 | 236.54 | 242.72 | 242.339 a 
135 122.747 122.49 | 125.17 | 124.532 = 
42 56.846 58.12 58.20 | 57.655 
17 24.315 25.68 24.76 | 24.417 a 
il 9.731 10.70 9.75 9.567 ag 
ll 5.457 4.47 5.42 5.267 <a 
x? with P| 20.294 34.52 25.52 | 25.939 eS 
v df. 8.d.f. 8.d.f 
= 
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K;/K, = 1 + 3(k — 1I)m+ (k — 1I)(k - 2)m’* — 3(K2/n) — (Ki/n’) 
and 
K,/K, = 1+ 7(k — 1I)m + 6(k — 1)(k — 2)m’ 
+ (k — 1)\(k — 2)(k — 3)m*® — 4(K;/n) — 6(K,K,/n’) 
— (Ki/n’) — 3(K2/nK,). 
The asymptotic efficiency in the restricted sense is 
E, = 1/D D.,,, . (28) 
The determinant of the covariance matrix up to order N~* for the 
sample zero frequency estimates of m and p is 
D(m, p) = m’as{a,(1 — Po) — PoK,}/N’azn’K,P, , (29) 
where a = 1 — p + p(l — m)*‘, a, = K./K, , Po = at and 
a, = 1 — (1 — m)* — km(1 — m)*". The asymptotic efficiency with 
respect to the method of moments is 
E, = D/D(m, p). (30) 
It may be noted that E, = Ey’ whenn = 1. 


TABLE 2 


Asymptotic Erriciency RELATIVE TO THE METHOD oF MoMENTS 
WHEN n = 1andDk = —-1. 


0. 
0. 
0. 
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INTRODUCTION 


In an earlier paper [1] we challenged the prevalent view that in 
comparisons of pain relieving drugs it is always desirable to have ‘“‘each 
patient act as his own control,” i.e., to test more than one drug on each 
patient and to estimate treatment contrasts from within-patient differ- 
ences. We sought data of other investigators to compare with our own, 
but we were unable to find references which gave sufficient detail to 
permit investigation of the merits of alternative designs and analyses. 
It. seemed desirable, therefore, to put our own data on record so as to 
facilitate discussion of the issues raised. This we did, and the response 
has been gratifying. 

This paper presents and discusses some of the suggestions and 
comments of those who wrote to us. We are particularly indebted to 
I'rederick Mosteller, Department of Statistics, Harvard University, and 
to Robert Curnow, A. R. C. Unit of Statistics, University of Aberdeen, 

Aberdeen, Scotland, who, in addition to making numerous thoughtful 
“comments, performed their own analyses of our data. 


A SIMPLE EXTENSION OF THE MODEL 


As stated in our earlier paper, this study of relief of post-operative 
pain was designed as an incomplete block experiment. Fach patient 
was given a dose of a test drug, and the dose was repeated when the 
patient again complained of severe pain. The number of hours of 
“greater than 50 percent pain relief’’ cumulated over both intervals 
was the measure of drug efficacy. A second drug was tested on each 
patient in the same way, starting at the time when the patient again 
complained of severe pain. Three drugs, two levels of a new drug and 
one level of Demerol, were under study. 
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It was found that the total hours of relief achieved when a given 
drug was the second administered was larger on the average than when 
that drug was administered first. This is consistent with the presump- 
tion that the pain decreases with time after operation. However, it 
invalidates a straightforward least squares analysis based on the model 


hours of relief = grand mean + patient effect 
+ drug effect + random error. 


We considered the possibility of extending the analysis to include a 
simple time-period effect, representing the average differences in relief 
between the second and first periods, but for two: reasons we did not 
pursue it. Firstly, and most importantly, it seemed to us that such a 
model would be unlikely to be an accurate reflection of the true situation. 
For example, if the duration of drug effect increases with time after 
operation, the period effect should be greater when the second drug is 
given after a potent drug than when it is given after a weak drug. 
Secondly, our object in the first paper was to compare simple alter- 
natives for the design and analysis of studies to evaluate analgesics, 
and a procedure requiring adjustment for time period as well as patient 
effects did not seem attractive as a procedure to recommend to 
clinicians for routine screening of new drugs. (In this connection one 
should remember that it is almost always necessary to discontinue study 
of some patients, often for reasons unrelated to drug response. Thus, 
balance in design is hardly ever achieved.) 

Curnow extended the least squares analysis of our data to take 
account of such a time-period effect. The model becomes 


hours of relief = grand mean + patient effect + drug effect 
+ time period effect + random error 


and the corresponding analysis is given in Table I. He stated that the 
differences due to the drug which follows a potent one being administered 
rather later than one which follows a weak one were small and had little 
effect on the estimates of drug differences. Thus the internal evidence 
of the experiment itself does not seem to confirm the fears expressed 
above, and the least squares analysis that includes a time-period effect 
appears to yield a reasonable description of the data. In particular, 
in both the inter- and intra-block analyses, the between-group residual 
mean square is now not larger than the corresponding within-groups 
mean square. 

‘Curnow went on to point out that the inter-block error is not much 
larger than the intra-block error, and that one might on this ground 
justify an analysis ignoring patient differences, as shown in Table II. 
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TABLE I 
Least Squares IntTRA-BLocKk ANALYSIS 


Source DF. 


Blocks 
Periods and Drugs 
Periods ignoring drugs 
Drugs adjusted for periods 
Between Groups Residual 
Within Groups Residual 


Total 


ComMPARISONS 


Mean Difference 


3.421 
2.037 
1.384 


TABLE II 
Least Squares ANALYsIS IGNORING PaTieENT DIFFERENCES 


Source S.S. 


Periods and Drugs 3 434.368 
Periods ignoring drugs 261.588 
Drugs adjusted for periods 172.780 

Remainder 82 982.888 


Total 85 1417 .256 


CompaRIsONs 


Mean Difference 


3.548 
1.861 
1.687 


578 
SS. MS. 
42 607.256 
3 378.015 
1 261.628 
2 116.387 58.194 
3 3.051 1.017 
37 428 .945 11.593 
| 85 1417 .256 
Comparison PS Variance t 
D vs. Ti 1.178 3.15 
D vs. T; 1.061 1.98 
T;3 vs. T; 1.025 1.35 
MS. 
86.390 
11.986 
Comparison Variance t 
D vs. Ti 0.874 3.80 
3 D vs.T; 0.831 2.04 
T; vs. T; 0.816 1.87 
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This observation suggests that in future experiments conducted under 
similar conditions we should not arrange matters to eliminate patient 
effects, if it will cost us much to do so. 

Our own calculations and scatter diagrams confirm Curnow’s re- 
sults. However, in view of our first argument against the use of a 
simple time-period effect we found them surprising and made some 
further investigation. We take note first that the fact that a model 
“fits”, in the sense that certain internal checks are satisfied, does not 
insure its correctness. If, for example, second-dose relief were more 
variable than first-dose relief, an examination of residuals would not in 
this case establish that fact; the estimated first-dose residual for a 
given patient has precisely the same magnitude as the estimated second- 
dose residual for that patient when we use the present model. The 
second-dose relief scores are, in fact, more variable than are those for 
first doses. This may be observed by examination of Table 1b in [1] 
and it is reflected in the fact that the within patient residual variance 
(Table 3a), which excludes the component of variability due to patient 
differences, is actually larger than the between patient variance for first 
dose relief (Table 2). 

If the facts are in accord with the assumption that the increment 
in relief with increasing time after administration is proportional to time 
after operation, the effect of the model having only a simple time-period 
effect would be to assign a part of the drug difference to the time-period 
effect and thus to reduce the measured differences between drugs, as 
compared to an analysis for the first period only. In fact, the difference 
between the most and least potent drugs (D and 7) estimated by use 
of the extended model is less by about 15 percent than the estimate 
obtained from analysis of the first period only. 

Granted, then, that the model with a simple time-period effect 
fits the data insofar as intra-block checks are concerned, what may we 
expect to gain or lose if we go ahead and use it? We do have some 
evidence which tends to contradict the assumptions of this model, but 
the effects do not seem large. We have some slight evidence that we 
may lose discriminating power by misinterpreting drug differences as 
time-period effects, but, as judged from these data, the reduction in 
variance may well compensate. Insofar as testing (i.e., deciding which 
drug is better) is concerned, or estimating potency in an experiment 
with a full range of standards, the analysis based on this model should 
give essentially valid results. 

In any event, apart from the question of the best analysis for these 
data, the smallness of the between patients component in the extended 
analysis suggests that we would have achieved at least equal precision 
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for drug comparisons if we had kept cach patient on a single drug 
throughout the experimental period. Had we done so, we would have 
had, in addition to possible benefits in precision, the comfort of closely 
imitating the ordinary clinical situation, and the advantage of an un- 
biased estimate of a clearly interpretable measure of effectiveness along 
with an unbiased estimate of its variance. We might also have had the 
cooperation of several more surgeons, some of whom, understandably, 
object to studies which require administration of several coded drugs 
to each patient. 


A REGRESSION MODEL 


In furtherance of the point of view that the cost of statistical work 
is generally small compared to the cost of gathering clinical data, 
Mosteller suggested that the data be analyzed in accordance with a 
model which might give a fairly realistic appraisal of the time-trend 
effect. In particular, he suggested that we not cumulate the two doses 
of each drug but analyze the data for single doses according to the 
following model. 


Yima + + Ymtima + Eimd 3 
where 


a; effect of subject, 
effect of medication m when given at time zero, 
= regression of effect of medication m on time of administration, 
time at which individual 7 receives the d-th dose of medica- 
tion m, 
€,ma = random error. 


This analysis was carried out, using the times of drug administration 
(not shown, but available upon request) in addition to the data already 
presented. It was found that the regression coefficients y,, were quite 
close to one another and, if we restrict the model to the case of equal 
Ym, the least squares estimate of this common value is 0.20, correspond- 
ing to an additional hour of drug relief for every five hour interval 
between operation and drug administration. Using this restricted model 
and doubling the estimated effects to make these results comparable 
to those for our other analyses, we have the estimates of treatment 
differences shown in Table III. 

Assuming for now the correctness of the regression model, we see 
that these estimates are quite close to those obtained in the analysis 
using only the first drug for each patient—Table 2 of [1]. However, 
the variances are reduced by about one-third, so that the regression 
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TABLE III 
CoMPARISONS FROM REGRESSION MopDEL ANALYSIS 


Comparison Mean Difference Variance t 


D vs. T; 4.066 0.824 4.48 
D vs. T3 2.530 0.732 2.96 
T3 vs. T; 1.534 0.708 1.82 


analysis, using all the data, appears to have about 50 percent greater 
precision than the “first drug only’”’ analysis, which uses only half 
the data. 

_ Comparing the apparent precisions provided by the above analyses, 
we see that Curnow’s intra-block analysis, adjusting for patients and 
time periods (Table I) appears comparable in precision to the “first 
drug only”’ analysis. Curnow’s analysis ignoring patient effects (Table 
II) appears almost equal in precision to the analysis based on the 
regression model. 


OTHER COMMENTS AND ERRATA 


A third correspondent, Irwin Bross, Roswell Park Memorial Insti- 
tute, Roswell Park, New York, raised several points, in part overlapping 
those above. In addition he pointed out that the results of the “‘least 
squares”’ analysis (Table 4a in [1]) differ appreciably from the analysis 
based on simple linear combinations of intra-individual contrasts. 
Since the design is nearly balanced, closer agreement might be expected, 
even though our “least squares” analysis fails to allow for time-period 
effects. In fact, contrary to the suggestion in [1], the “linear combina- 
tions” analysis is not really comparable to the ‘“‘least squares’ analysis. 
The reason is that in estimating, say, D — T, , no account was taken 
of the information about this difference given by contrasting the direct 
estimate of D — T; with the estimate of T, — T; . If we take an 
average of the directly observed difference, D — T, , with the contrast 
obtained by combining D — T; with T, — T; , weighting inversely as 
the estimated variances, and proceed similarly for the other compari- 
sions, we get Table IV, which agrees much more closely with the “least 
squares’’ analysis (Table 4a in [1]). As might be expected, the “ex- 
tended least squares’ analysis which does allow for time-period effects 
(Table I) gives results in excellent agreement with Table IV. 

This revision leads in turn to a more precise combined intra- and 
inter-block analysis in place of that shown in Table 3c in [1]. The 
result (not given) is quite close to the “least squares’ combined analysis. 
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TABLE IV 
DrucG ComMPARISONS FROM ANALYSIS OF WITHIN PATIENT CONTRASTS 
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Comparison Mean Difference Variance t 
D vs. T; 3.449 1.206 3.18 
D vs. T; 2.051 1.073 1.98 
Ts vs. T; 1.398 1.031 1.38 


w 


TABLE 4c 


1 
0.8518 


(revised ) 


w’ 


CoMBINED INTRA- AND INTER-BLOCK ANALYSIS FOR DruG DIFFERENCES 


~ 3.720 
Calculation of Average Mean Difference: Drug D vs. Drug 7, 


Both Curnow and Bross pointed out that erroneous entries are 
shown in Table 4c in [1]. This table should be replaced with Table 4c 
(revised) as shown. 


0.2688 


(w)(3.150) + (w’)(3.900) _ 


w+w 


1 1 


w+w 0.8518 + 0.2688 


DrucG ComPARISONS 


Calculation of Variance of Difference 


= 0.892. 


Mean Difference 


.330 


3 
1 
1 


ooo 
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DISCUSSION 


The point.which we wished to emphasize in [1] was that on account 
of the increase in duration of relief with time of administration, the use 
of simple intra-patient contrasts—‘‘each patient his own control’’—may 
not be optimal, and that an experiment in which each patient received 
only a single drug might be preferable. 
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The analyses proposed by Curnow and Mosteller have demonstrated 
that for this data a model taking account explicitly of the effect of time 
of drug administration will yield more efficient estimates than does the 
one-way analysis of the first period only. With respect to the analysis 
of the present data, both the Curnow and Mosteller models appear to 
fit quite well, and the estimates of the drug effects produced by their 
analyses are more precise than those derived from the analysis for the 
first period only. 

It is not clear, however, whether the above finding should be con- 
strued as evidence supporting the need for more complex models and 
analyses than are currently in vogue, or whether instead it is evidence 
supporting our original viewpoint. If, as the evidence of Curnow’s 
analysis suggests, patient differences are negligible, we would have no 
reason at all to use each patient as his own control. Were each patient 
given only one drug, the time-period effect of Curnow’s model would be 
eliminated from treatment contrasts and the one-way analysis—using 
now the whole of the data—should be fully efficient. 

Even if patient differences are not negligible, the advantages of 
simplicity, boi of design and analysis, may outweigh a small gain in 
efficiency which could be achieved with the use of a more complex 
design and analysis. With our data, for example, the estimates of 
greatest apparent precision were those provided by the regression 
model, and the gain in precision compared to the analysis of the first 
period only was 50 percent. Such a gain is, of course, quite worth 
having, but it is not overwhelming, and it is easy to believe that had 
our design assigned one drug to each patient throughout the period of 
study, the one-way analysis might have equalled the present regression 
analysis in precision. 

In conclusion, we must: again emphasize that we do not claim that 
our findings apply to all kinds of pain studies. There may be many 
situations for which within-patient contrasts are far superior to between- 
patient contrasts. What we do claim is that no principal such as 
“each patient his own control’ is entitled to the status of dogma. 
In some situations, at least, the simpler methods are better. 
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FITTING A GEOMETRIC PROGRESSION 
TO FREQUENCIES 


E. J. WILLIAMS 
Division of Mathematical Statistics, 
C. 8.1. R.O., Canberra, Australia 


SUMMARY 


This paper discusses the interpretation of frequency data when the 
series of observed frequencies is assumed to arise from a population in 
which the expected frequencies form a geometric progression. Such 
situations occur in the study of steadily increasing insect populations 
and similar phenomena, where the common ratio of frequencies is related 
to the rate of growth of the population. 

The estimation of the common ratio, and tests for the significance of 
departure from geometric trend, are discussed. Asymptotic formulae 
for the tail probabilities in large samples are determined. 

The methods and tests of significance are illustrated by siiiliaatias 
to some experimental data. 


I. INTRODUCTION 


Observations are often recorded by classifying them into several 
classes and counting the frequencies of occurrence in each class. In 
general, the expected frequencies will be partly specified by theoretical 
considerations, but will often also depend og one or more unknown 
parameters. 

The interpretation of such data then involves, firstly, testing its 
concordance with the assumed form of expected frequencies, and 
secondly, estimating the unknown parameters. Thus, for instance, in a 
steadily increasing population of organisms, the number of individuals 
expected in successive age-groups will be in geometric progression, the 
common ratio between the expected numbers representing not only the 
growth of the population but also the effects of mortality, migration 
and other factors. In such a study, one of the objects would be to test 
the concordance with the assumed geometric progression, and the other 
would be to estimate the common ratio accurately. 

In general problems of this kind, the parameters may be estimated, 
and the accuracy of the estimates assessed, by the method of maximum 
likelihood. To test the fit of the model to the data, the x’ test is generally 
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appropriate. If hypothetical values of the parameters have been 
specified, as occurs in many practical instances, their concordance with 
the estimates derived from the data ‘can also be tested by means of x’. 

In the rarticular case, with which this paper deals, of expected 
frequencies in geometric progression, there exist sufficient statistics for 
the two parameters involved, representing the general level of the 
frequencies and their geometric trend. This facilitates estimation and 
significance testing, although, since the equations of estimation are in 
general non-linear, iterative methods are required in the arithmetical 
work of calculating estimates. 

The need to fit frequencies in this way arose from a method, devised 
by Dr. R. D. Hughes, Division of Entomology, C.S.1.R.0., of using 
the age distribution in an aphid population to study its rate of growth 
in the field. The device of using the proportions of individuals in the 
immature instars in a field population of aphids to give, firstly, the age 
distribution, and secondly, the rate of increase of the population, is 
the subject of a forthcoming paper by Hughes. In the present paper, 
some preliminary observations of instar distribution under controlled 
experimental conditions are discussed. 

Chapman and Robson [1960] have considered the age distribution in 
a stationary population subject to constant mortality. As this formula- 
tion leads to a geometric progression of expected frequencies, many of 
their results anticipate results given in this paper. However, the 
general objects and scope of the two papers are different. 

Aphid populations under uniform and favourable conditions increase 
at a constant rate, so that their expected numbers at equal intervals 
of time form a geometric progression. However, the numbers even of 
an initially small population quickly become too large to be counted 
accurately. An estimate of the rate of growth of the population may 
then be made from the age distribution. The immature aphid passes 
through four instars before reaching maturity. In a stationary popula- 
tion, the number expected in each instar is proportional to its duration. 
For some species of aphids, for instance A. craccivora, the average 
durations of the first three instars are probably equal under constant 
conditions, so that equal numbers will be expected in each instar in a 
stationary population. When the population is increasing at a constant 
rate, the numbers in each instar will approximate to a geometric pro- 
gression. 

If the common duration of each of the first three instars is c, and 
the growth-rate of the population is p, then the ratio of the number in 
each instar to that in the preceding one is approximately 
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In this result, mortality has been neglected, since it has been shown 
that under favourable conditions mortality in the immature stages is 
negligible. 

By taking a random sample from the aphid population and de- 
termining the number in each instar we can first check the validity of 
the assumption of uniform growth by testing whether the numbers are 
consistent with a geometric progression. Having satisfied ourselves 
on this point, we can then estimate the common ratio of the numbers. 
If c, the duration of each instar, is known, the growth-rate can then 
be determined. 

If the durations of the different instars are different, the above 
method will be modified, and the estimation of the growth-rate is a 
little more complicated; in principle, however, such data will still 
provide information from which the growth-rate can be determined. 


II. THE BASIC DISTRIBUTION 


We consider a sample of N individuals classified into k classes, the 
observed frequency in class 7 being n; , with expected value 


E(n,) = @=0,1,---,k—1). 


If N is assumed to have a Poisson distribution, so has each of the n, ; 
the joint probability density is then 


—Ado (u) 


P(m , Il Mo !m ! +++ my!’ 


(1) 
where 
and 


From the form of the density it is apparent that N and T are a pair 
of sufficient statistics for the parameters \ and yw. Thus, questions of 
estimation may be referred to the joint distribution of N and T. 

By summing the probability (1) over values of the n; leading to the 
given values of N and 7, we find the joint distribution of N and T as 


P(N, T) (2), 


where C(N, T) is a numerical coefficient. 
The number WN has a Poisson distribution with mean A¢o(u); 


P(N) = 


1 
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Hence for given N, the conditional density of 7, which is independent 
of X, is 


P(T | N) = C(N, 


This is a special case of the distribution discussed by Noack [1950] and 
described by him as a power-series distribution. Our [¢(u)]” takes 
the place of his f(z). From this expression for the density we see that 
C(N, T) is the coefficient of u” in ¢o(u)”. This fact enables the numerical 
coefficient to be determined directly. 

This conditional density of 7 provides the basis for estimates and 
tests of significance about pz. 


Ill. THE COMBINATORIAL FACTOR 


The coefficient C(N, T) is seen to be a generalization of the binomial 
coefficient (for the positive binomial, when k = 2; for the negative 
binomial, when k = ). It is the number of ways of allocating T like 
objects among N different cells, none of which may contain more than 
k — 1 objects (see Riordan [1958], page 104, where some recurrence 
relations and moment formulae are also given). 

The coefficients are generated by the function 


On equating coefficients of u” we find 


own = 


The series has 1 + [7'/k] terms. This expansion in terms of binomial 
coefficients is useful, especially if k is large, since the first few terms 
then give a close approximation. 
Because of symmetry, C(N, T) = C[N, (k — 1)N — T] 80 that 
results for T > $(k — 1)N can be found from those for T' < $(k — 1)N. 
From the generating function or otherwise we may also deduce the 
recurrence relation 


C(N,T) = C(N,T —1)+ C(N —1,T) —C(N —1,T — &) 
which is useful for computation. 
Hitherto we have been considering relations among the C(N, T) 


for a fixed value of k. We now consider relations for different k, and 
shall indicate by a suffix the number of classes. 


+ 
4 


= 
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We may express the coefficients for k classes in terms of the co- 
efficients corresponding to the factors of k, by means of the following 
reduction formula. 

Ifk = fg, 


1- 
= 


Hence 
CAN, T) CAN, T) +> CN, ACAN, 1) 
+ CAN, T — 2) + 


a series of 1 + [7'/f] terms. The terms of these series are all positive, 
and generally smaller than those of the series (3) of products of binomial 
coefficients. They may have some advantages for computation, pro- 
vided the coefficients for f and g classes are known. 

By differentiating the generating function with respect to its param- 
eter we may prove recurrence formulae of the type 


N 

— — 1,T) + & — 3)CW 1,T — 1) 
+ — OCW — 1,T —& + 


When k = 3, many particularly simple recurrence formulae may be 
established. Recurrence relations when k = 3 are 


— 1)C(N — 1,T — 1) +3(N — 1I)C(N 2,T — 
1 


2T 


[(2N — T+ 1)C(N,T — 1) + 3NC(N — 1,T — 2)] 


= lew - 1,7) CW - 1,7 - 


the last being a particular case of (4). As k increases the recurrence 
formulae become more complicated. 


IV. MOMENTS OF THE CONDITIONAL DISTRIBUTION OF T 


The moments and cumulants of the conditional distribution of T 
are of general interest, and will also be of use when the determination 
of probabilities in the tails of the distribution is being considered. 


| 
T(2N — T) 
| 
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Since the cumulants for a sample of N are simply N times those for a 
sample of 1, we shall consider a sample of 1, for which 7 equals X, the 
number of the class (0, 1, --- , & — 1) in which the observation falls. 


Clearly, if 
( d k-1 
o-(u) = do(u) = 
then the rth moment about zero is ¢,(u)/¢o(u). This result is equivalent 


to the results of Noack [1950], though expressed in a slightly different 
manner. 


However, unless k is small, these expressions. are not convenient 
for the computation of the central moments and cumulants, which are 
more simply found directly, 


The moment-generating function of X is 


E@’*) = = do(ue’)/bo(u). 


Thus the cumulant-generating function is 
K(s) = log do(ue") — log do(u) 
= —log (1 — ue’) + log (1 — ue) + log (1 — u) — log (1 — »’) 
r rk 


From this expansion we derive the particular results 


ky" 
1 1 — 1 ’ 
k? k 


= 


(1 — (1 — 


The form of the cumulant-generating function shows that, if u, is a 
geometric variable with parameter u*, then X + ku, = u,. It also 
follows that, if u,(N) is a negative binomial variable with parameter y* 
and index N, then T + ku,(N) = u,(N). 

When pu = 1, the distribution is the uniform discrete distribution, 
and the cumulants are expressible in terms of Bernoulli's numbers. 
We then have 


k 2k 
2 3 a 
Kg = 
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= B, s* 


2 
where s/(e’ — 1) = 1+ /r! = e** in symbolic form. Hence 

= (k” — 1)B,/r; in particular, x, = (k — 1)/2, x, = (k° — 1)/12, 
fi. = 0, i, = — (k* aa 1)/120. 


Relation between Cumulants Corresponding to 4 and yw" 


When the series of frequencies is reversed, the ratio yu is replaced 
by its reciprocal. It therefore follows that, when yu is replaced by u™', 


the odd cumulants other than «x, are simply changed in sign, and the 
even cumulants are unaffected. We have 


=k —1—x«,(y), 
k(u') = (—1)"«(u) (r > 1). 


These results may also be readily verified from the form of the cumulant- 
generating function. 

Because of these facts, we need not consider values of » outside the 
range (0, 1). 


V. ESTIMATION OF THE RATIO OF FREQUENCIES 


In the conditional distribution, 7 is a sufficient statistic for the 
parameter yu, the common ratio of the expected ferquencies. Thus the 
estimation of u is straightforward. However, since the equation of 
estimation by the method of maximum likelihood is non-linear, iterative 
methods will be required to solve it. 

The logarithm of the conditional probability of 7’, apart from terms 
independent of the parameter, is L = T log » — N log ¢o(u), and its 
first derivative with respect to yu is 


(T/u) — (5) 


Equating the derivative to zero gives the maximum likelihood estimator 
of u. We indicate by an asterisk the maximum likelihood estimator and 
functions of it. Thus 


— = (6) 


Since the equation is linear in 7’, it is clear that the estimator is to be 
found simply by equating T to its expectation: 7 = Nx* , as may be 
verified from the results of the previous section. 
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The variance of the derivative (5) gives the information I » about pz 
in the sample. Being linear in 7, the derivative has variance 


This result is a special case of the results of Patil [1961], who investi- 
gates the estimation of the parameter of a generalized power-series 
distribution. In large samples, the reciprocal of J, approximates the 
variance of the estimate u*; that is, V(u*) = w’/Nxe. 

For purposes of calculating and tabulating the solutions of the 
maximum likelihood equation, we shall put 


T/(k — 1)N =», 


so that, for all k,O < v < 1. 
Then the estimating equation may be written 


7 = 1)v. (7) 


An alternative form is 


l—u 


= 1)(1 = v), 


ky * 


Thus, if u* is the root corresponding to v, then u*~’ is the root 
corresponding to 1 — v. In particular, ifv = 4, u* = 1. 

Since the roots corresponding to values of v exceeding } are the 
reciprocals of roots corresponding to values of v less than 4, we may 
henceforth confine attention to v < 3, 4 < 1. 

Equation (7) can be solved iteratively, once an approximate value 
of » has been chosen. If v is not too near 4, and k is not small, a first 
approximation is 


and a second approximation is 
= + kui(l — ws)’. 


The difference between the two sides of (7) when an approximate 
value of u is substituted represents —u/N times (5), the first derivative 


| 
or 
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of the likelihood. The adjustment to yu is given by the ratio of this 
difference to —uJ,/N. Then, approximately, 


The substitution and adjustment may be repeated as often as required 
to give the desired accuracy. 

If v is close to 3, an approximation alternative to (8) is to be pre- 
ferred. We put 


p=e’, v=}—w. 


Then equation (7) becomes 


(k —w 
1-—e”’ l1—e MG 
or, symbolically in terms of Bernoulli’s numbers, 
1 


We then find 


— 1) + — 1) + 06) = (k — 1)w, 


whence 
and 
12w 4 72w* 144 (k? + 11)w* 


5 + OW). ©) 


Solutions of the maximum likelihood equation for various values 
of k and v are given in Table 1. Once values of u are given, it is easy 
to compute J, . As we shall see in Section X, the solution not only 
gives a point estimate of the ratio 4, and an approximate standard 
error based on the information function J, , but also gives a means of 


determining tail probabilities, needed in making significance tests and 
setting confidence limits. 


VI. ESTIMATION OF 


In general, for the problems considered in this paper, the actual 
value of X, representing the size of the population (or rather, of the 
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TABLE 1 

Maximum EstiMaTE OF [WITH ARGUMENT v = T'/(k — 1)N] 

k 2 3 4 5 
0.00 0 0 0 0 
0.05 052632 092894 131334 167119 
0.10 111111 178395 238358 291010 
0.15 176471 261940 333333 392939 
0.20 250000 346500 422530 483323 
0.25 333333 434259 509668 567737 
0.30 428571 527202 597388 649654 
0.35 538462 627431 687922 731617 
0.40 666667 737405 783468 815815 
0.45 818182 860221 886484 904443 
0.50 1.000000 1.000000 1.000000 1.000000 


zero-class), is not of interest. However, for testing the significance of 
departure of the frequencies from a geometric progression, estimates of 
the expected frequencies in each class will sometimes be required. 
We shall then be interested in the simultaneous estimation of \ and yu, 
based on N and 7. We therefore consider this question, rather than 
the estimation of \ alone. ‘ 

The joint probability of N and 7, as given in (2), is 


C(N, T) 
N! 
so that the logarithm of the likelihood, considered as a function of the 
two parameters, and with constant factors ignored, is L = —)o(u) + 
N log\ + T log u. The derivatives are 


L, = [—Adi(u) + 


The estimates are given when these derivatives are equated to zero, 
or when N and T are equated to their expected values. For direct 
solution, X can be eliminated between the two equations, giving an 
equation for » identical with (6). « having been found, \ may then be 
found by substitution in (10). 

For simultaneous solution, and also to give a measure of the infor- 
mation about A and yu, we require the covariance matrix of the deriva- 
tives. Noting that N is a Poisson variate, with variance equal to its 
mean A¢o(u), that the covariance of N and T is \¢,(u) and that the 


4 
= —¢o(u) + (N 
2 
| 
. 


594 BIOMETRICS, DECEMBER 1961 


variance of T is \¢2(u), we find for the covariance matrix of the deriva- 
tives 
Loi /Au 
The inverse of I, giving approximate variances and covariance of 
the estimates in large samples, is 
— 


From trial values \, » of the parameters we find improved estimates 
hy means of the equations 


which give explicitly 


(11) 


= 


| 
— 


No. — |. 
— $1) 

Note that the trial value of \ does not appear in the equation for 
\*, since the estimating equations are linear in \*. The adjustments 
given by (11) may be repeated to give the accuracy required. 

The expected frequency in class 7 is \u'; the variance of the esti- 
mated frequency may be determined, using I~’, as 

2i 
dud. — 


The relative variance of the estimated frequency is therefore 


o2 — 2i¢, + 
— $1) 


VII. TESTS OF DEPARTURE FROM PROPORTIONALITY 


(6. — 2i¢, + 


Before any use is made of data to which proportional frequencies 
have been fitted, it is necessary to test whether the observations depart 
significantly from the estimated frequencies. 

Following Cochran [1954], we shall denote by X’ the statistic used 
for testing the discrepancy between observed and expected frequencies, 
and by x’ the random variable with the familiar distribution, to which 
the distribution of X* approximates. We consider here the adequacy 


LL. 
= No: — Td, 
— $1 
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of the x’-approximation in testing X’ for departure from expectation 
in the present problem. 

Alternatively, we may use the fact that, because N and T' are 
sufficient for \ and y, the probability of a sample, conditional on N and 
7, is independent of the parameter values. An exact conditional test 
can thus be made, provided there are several possible samples corre- 
sponding to the given values of N and 7. The probability of a sample 
of values no, , With = N and in, = T is 


— 1 T 
Likewise the joint probability of N and T is 


where C(N, T) is as defined in (2). Hence, the probability of the sample 
conditional on N and T is 


N! 


By enumerating all possible samples and their probabilities we can 
decide the significance of an observed sample. Often the set of samples 
deemed to show significant departure from proportionality will be taken 
as a set of the samples with the smallest probabilities; however, in 
many problems, other criteria may be considered more appropriate. 
The use of X* has the advantage that it orders the samples according to 
a quantitative measure of departure from proportionality, whereas the 
probability attaching to a particular sample may be small merely be- 
cause of the discontinuity of the distribution. A satisfactory arrange- 
ment would therefore be to order the samples by means of X’, but 
to use the exact probabilities. In many cases the exact cumulated 
probabilities will not differ much from those given by the x’-distribution. 
The questions of calculation of probabilities and of choice of the signifi- 
cant set have been thoroughly discussed by Fisher [1950]. 

As an example of the exact calculation, we consider the following 
data in four classes: 


Observed Expected 


30 
5 
2 
3 


no 27 
ny 9 
3 
: 
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It will be seen that, because the numbers are restricted to four classes, 
the possibilities are much fewer in number than in other examples, 
such as that of the testing of deviations from a Poisson distribution, 
as discussed by Fisher. For this set, N = 40, T = 18, and it is readily 
verified that u* = 4,\* = 27. X’ is easily calculated as 


2 2 2 2 
N= 64M; 


however, since the expected values in two of the classes are small, 
the tabulated distribution of x’ is not likely to be a good approximation 
to the distribution of X’ for this sample. For this reason, as an example 
of method, the exact probabilities have been determined for the 37 
possible samples having values N = 40, 7 = 18. The probabilities 
and the values of X’ are given in Table 2. 

We have C(40, 18) = 220,495,831 x 10°. 

From the Table it is seen that the sample in question just attains 
significance at the 5 percent level; this is so whether the samples are 
ordered according to their probability or by X*. The tabulated distri- 
bution would give a probability, for x* with two degrees of freedom, 


e = — 0.039866. 


This shows that, for this example, the tabulated distribution under- 
estimates the significance probability only slightly. On the other hand, 
for the sample 


for which the probability accumulated according to X’ is 0.008190, 
= 8.593, 
the probability of exceeding which is 
0.013619. 
The probability here is overestimated somewhat. __ 

As a matter of interest, the mean and variance of for the distribution 
generated by this set are E(X’) = 1.9956, V(X’) = 3.4223, compared 
with the values of 2 and 4 respectively for the x’-distribution. The 
low variance will usually lead to overestimation of probabilities at the 
: tails of the distribution and underestimation of significance, when 


referred to x’, though this effect will sometimes be masked by local 
irregularities of the distribution. 


VIII. MOMENTS OF THE NON-PARAMETRIC DISTRIBUTION 


We now show how the moments of the exact distribution of X? 
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TABLE 2. 
THE FREQUENCIES IN Four CLASSES FOR N = 40, T = 18 
Cumulative 
no ny Ns Probability probability 27 x? 
34 0 0 6 .000000 -000000 1048 
31 0 9 0 -000001 -000001 610 
33 0 3 4 .000003 . .000004 522 
33 1 1 5 .000004 .000008 696 
32 0 6 2 .000010 .000018 376 
32 3 0 5 .000019 .000037 646 
31 1 7 1 .000089 .000126 352 
32 1 4 3 .000098 .000224 334 
32 2 2 4 -000146 .000370 424 
30 2 8 0 .000173 .000543 408 
22 18 0 0 .000514 .001058 376 
31 4 1 4 .000781 -001839 370 
30 6 0 4 .000807 .002646 360 
31 2 5 2 .000938 .003584 226 
31 3 3 3 .002083 .005667 232 
30 3 6 1 .003229 .008896 198 
29 4 7 0 .003460 012356 250 
24 15 0 1 .004561 .016917 198 
28 9 0 3 .005574 .022491 190 
23 16 1 0 .006841 .029332 226 
26 12 0 2 .009578 .038910 136 
30 5 2 3 .009688 .048598 174 
30 4 4 2 012110 .060708 120 
29 7 1 3 .013840 074547 160 
28 6 6 0 .023412 097959 136 
29 5 5 1 .029063 127022 88 
24 14 2 0 .034206 . 161228 120 
25 13 1 1 .038311 199539 88 
27 10 1 2 .046824 . 246363 66 
29 6 3 2 .048439 . 294802 58 
27 8 5 0 070236 .365038 66 
28 8 2 2 .075253 .440291 40 
25 12 3 0 .083006 . 523297 58 
28 7 4 1 . 100337 .623634 22 
26 10 4 0 . 105354 . 728988 40 
26 11 2 1 .114932 .843920 22 
27 9 3 1 . 156080 1.000000 0 


may be determined. This distribution is independent of the population 

parameters, and depends only on the values of N and T. 
We consider first the conditional moments of the n,; . 

pression (12) we see that 


From ex- 


‘ 
i } ; 
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where the sum is taken over all sets of the n; such that > n, = N and 
> in; = T. Now the ig value of n; is 


In the sum, the total of the arguments, with n; replaced by n; — 1, 
is N — 1, and the weighted sum of the arguments is 7’ — j. It follows 
that E(n;) = NC(N — 1, T — j)/C(N, T). By a similar method it 
follows that, if }> a; = A, and )> ia; = B, then the general factorial 
moment is 


Thus in principle the joint factorial moments of the n; can be de- 
termined provided we can regard the C(N, 7) as known. Some dis- 
cussion of the asymptotic representation of C(N, T) is given in Section 
X, but more needs to be known about these coefficients in general. 

For the test of departure from expected frequencies we have 


2 
x* = N. 
Now \* and y* are functions of N and T only, so are fixed for the con- 
ditional distribution. Hence X* may be simply regarded as a weighted 
sum of squares of the n,. In particular, since 


NW — — 2,T — 2%) + NC(N — 1,T — 1) 
= T) 
we have 
E(X’) = N 


A*C(N, 


(N — IC(N — 2,7 — 21) + C(N — 1,T — 2) 


higher moments may be found similarly. 


IX. TEST FOR TREND IN FREQUENCIES 


Of particular importance is the test for the reality of any apparent 
trend in the frequencies—that is, whether the estimated value of u 
differs significantly from unity. When N is large, the test may be 


| 
| 
\ 
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made by comparing the difference 1 — y* with its standard error, or, 
what is equivalent, testing $(k — 6d — T against its standard error. 
When uz = 1, the variance of u* reduces to 12/ (k? — 1)N, so that 
the test statistic, distributed nent as x’ with 1 degree of 
freedom, is 
(k? — I)N(1 — w*)*/12. (13) 
More directly, the variance of T is 
— 1)NJ/12, 
so that an alternative test statistic is 
12(4(k — —TY/(k DN. (14) 


The statistics (13) and (14) are equivalent in large samples since, 
when u* is near to unity, we have 


12(4(k — 1)N — T) 


=1— DN + O(N"), 


as can be seen from (9). 

The exact significance probability is given by the probability of a 
value of T less than or equal to that observed, which is easily computed 
directly when N is not large. This probability should be doubled, 
since both tails of the distribution are relevant to this test. With 
u = 1, this probability is simply 


C(N, = S(N, T)/k*. 


Now the generating function of S(N, T) is 


— w) = (1 — — 
From this representation we readily deduce that 


sw) = ("5 7)- 


which provides the most convenient means of computation of isolated 
values of the sum. On putting 7 = (k — 1)N in (15) we have the 
interesting corollary 


it 
1 r=0 

(15) 

ae 

at 

ne 
a: 
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Values of 7 less than 4(k — 1)N required for significance at the 5 
and 1 percent levels of probability have been determined for various 
values of k and N, and are presented in Table 3. For values of N beyond 
the range of the Table, the large-sample x*-test may be used. 


TABLE 3 


Exact ror Exisrence or TREND 
(Le. or Nutt. Hrrorsesis = 1) 


Panconr Pours or T (k — 1)N — T] 


N 3 4 5 
5 1 2 3 
6 1 3 4 

7 2 4 6 
8 3 5 7 
9 3 6 9 
10 4 7 10 
il 5 8 12 
12 6 9 13 
13 6 1 15 
14 7 12 17 
15 8 13 18 


1 Percent Point or T [orn (k — 1)N — T] 


= 
w 
ao 


Bebo 
OO 


= 
4 
at 
Ty 
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X. EVALUATION OF TAIL PROBABILITIES 


When N is large, the calculation of the probabilities of the tails 
of the distribution is troublesome. Blackwell and Hodges [1959] have 
given a method for expeditiously finding approximate tail probabilities. 
The method depends on a transformation of the distribution to a 
new distribution for which probabilities in the neighbourhood of the 
mean are to be determined. Daniels [1954] has given equivalent re- 
sults, expressed explicitly for continuous probability densities, and 
Good [1957] has given similar results to the term of order N~* of the 
leading term. 

We briefly outline Blackwell and Hodges’ results applied to the 
present distribution. We make use of the fact that the distribution 
of T is the N-fold convolution of the distribution of X as defined in 
Section IV. 

The moment-generating function of X — a, where a is any constant, 
is M(s) = Efe**~]. 

Let s* be the unique value of s that minimizes M(s), and denote 
the minimum of M(s) by m(a); note that s* is a function of a. Then 
we transform to a new variable Y, for which 


P(X 
m(a) 


P(Y = 2) = (16) 
It is readily seen that Y has in fact a probability density with the same 
range as X. Also, by differentiating the numerator of (16) with respect 
to s* and summing, we find that E(Y) = a. 

Before continuing with Blackwell and Hodges’ method, we establish 
the interesting result that, for any distribution for which the sum of 
the values of X is a sufficient statistic for the parameter yu, the new 
variate Y has a distribution of the same form, but with parameter u*, 
the maximum likelihood estimate corresponding to the value X = a. 
This result is probably well known, though we have never seen it 
discussed. Thus we see that, for distributions of this commonly occurr- 
ing class, the solution of the maximum likelihood equation is important 
not only for providing a point estimate, but also for evaluating the tail 
probabilities corresponding to the particular observed value. 

If the sum of the values of X is sufficient for », the density of X 
may be written as 


h(X)e""/f(u), (17) 


where 


= 
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Then 
M(s) = /f(u), 


-o that, for any value of s whatever, a transformed density is given by 


‘This density is clearly of the same form as that of X. 


Now the maximum likelihood equation, corresponding to an observed 
value a, turns out to be 


th(x)e” / >> 0, 


where we have differentiated (17) with respect to g(u) rather than yp 
itself. This is in fact F(X | u» = u*) = a. Now from (16) we derived 
that L(Y) = a. It follows therefore that the parameter of the distribu- 
tion of Y is p*. 


We can also express m(a) in terms of the maximum likeli- 
hood estimator. On putting « = a in (16), we see that m(a) = 
P(X = a)/P(Y = a), or Pa | w)/P(a | 

lor the particular distribution we have been considering, therefore 


mo) = (2) 


where we indicate by an asterisk a function of yu*. 


Blackwell and Hodges, by means of the transformation (16), and 
taking Na = T, show that 


P(X, + Xy =T) = [m(a)]*P(Y, + + + Y,=T), 


and thence, that each side is equal to 


T N 2 


(18) 


where the «x* are the cumulants of the distribution of Y; that is, fune- 
tions of u*. . The explicit terms in (18) give a satisfactory approximation 
to P(T | N), provided N is large. 

We note that the series in (18) depends only on y* (i.e., on 7), and 


not on w. Thus what we have obtained is an asymptotic expansion of 
C(N, T). In fact, 


A 
z 
z z 
or 


FITTING A GEOMETRIC PROGRESSION 603 


C(N, T) = E + (4, + ow} 
Thus, C(N, T) may be estimated by substituting the actual values of 
the cumulants, from the formulae given in Section IV. 
As a check on the asymptotic formula, we estimate C(N, T) for — 
the values k = 4, N = 15, T = 11. We deduce u* = 3, ¢ = 15/8, 


194/15” = 0.862222, 
2842/15° = 0.842074, 
1434/15‘ = 0.028326. - 


* 
2 


K. 


* 
3 
* 
4 


K 


Hence 


1 


The true value is 
_ (18)(17) 
from which the estimate by the asymptotic formula is in error by 130 
per million. 
For significance testing, the cumulative tail probabilities are more 


important than the individual terms. These cumulative probabilities 
are also given by Blackwell and Hodges. They give 


_PT=)[, 2 +8042 4 ] 
E NN OWN) 


where, for the distribution considered here, z = u/u*. We note that, 
as is to be expected, the bracketed series in the cumulative probability 
depends on yu. The approximation is valid only for u* > yu, and is 
most effective when z is small—that is, when u* is large and the extreme 
tails of the distribution are being considered. 

When we are determining the probability in the lower tail of the 
distribution, that is P(T < t), so that T < $(k — 1)N, and yp* < gp, 
we simply replace z by z' throughout the formula, and change the 
sign of x«*. 


j 
C(15, 
Bet 
= 
{ 2 
ait 
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As a test of the adequacy of this approximation, we apply it to the 
example for which the individual term was calculated above, namely 
k=4,N = 15,7 = 11. 

Inserting the values found above, we have z = 2u, 


P(T = 11) 


(, 1) 4 4 2) 
1 ( 2) + 
On 194° (1 
15* 
We consider the adequacy of the approximation when » = 1, for which 
we can readily determine the exact probability. We have 
P(T = 11) = 2,784,600/4"°, 
so the leading term in the approximate probability is 
2P(T = 11) = 5,569,200/4"° = 0.005187. 


For the first adjustment we find the factor 1 — 0.078223 giving 
5,133,560/4"° = 0.004781. The exact value of the denominator is 


giving the exact probability 0.004893. 

We see that the first adjustment gives an approximation adequate 
for the purpose of assessing the significance probability. Since both 
tails of the distribution are relevant to any test of significance, the 
significance probability is 0.009786; 7’ = 11 is thus the 1 percent point 
of the distribution. _ 

By determining the tail probabilities corresponding to any given 
value of » we can in principle set a confidence range for » — the set of 
all » which are not discordant with the observed values N and 7’. 


PT < 11) = 


+ O(N~%) |. 


XI. APPLICATION TO EXPERIMENTAL DATA 


The data in Table 4 show numbers in each of the first three instars 
of immature cowpea aphids (A. craccivora), from samples drawn on six 
different occasions from the same population. Since the experimental 
conditions were uniform, and the average duration of each instar is the 
same (about 42 hours at 20°C), the expected numbers are in geometric 
progression. The Table also gives an overall estimate of the common 
ratio, 0.529820. In order to determine the growth-rate of the population, 


ud 
| 
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TABLE 4 
AcE-DIsTRIBUTION OF IMMATURE ApuHips (A. craccivora) 
Instar 

Sample I II Il N T x 
1 181 92 51 324 194 179.730 0.526015 | 0.11 

2 148 78 42 268 162 147.737 0.531518 | 0.01 

3 130 54 43 227 140 123.456 0.543412 4.07 

4 70 37 21 128 79 69.580 0.543848 | 0.03 

5 88 41 20 149 81 87.710 0.474049 0.13 

6 85 42 28 155 98 82.857. 0.558629 0.63 
Total 702 | 344 | 205 | 1251 754 690.958 0.529820 2.14 


we equate this ratio to e”, c being the average duration of instar and 
the growth-rate. 
We then find 


cp = —log (0.529820) = 0.6352 
so that 
p = 0.6352/42 per hour = 0.01512 per hour = 0.363 per day. 

The population grows at the rate of 100(e”'*** — 1) = 44 percent per day. 

Table 4 also gives estimates of \ and u for each sample, and X’ 
with 1 degree of freedom measuring departure from the common ratio. 
Only the largest X*-value, for sample 3, attains the 5 percent level of 
significance, so that there is no evidence for departure from hypothesis. 


The variance of the estimate of » from a sample of N is u’/V(T) 
which reduces when k = 3 to 


+ 
N(1 + 4u + w’) 


Using the overall estimate ut = 0.529820 we find the variance of esti- 
mate to be 


V(u* | N) = 0.510813/N. 


The consistency of the estimates of » from the different samples 
may be tested by X’. 
For consistency, 


X? = (0.526015’ 324 + 0.531518 268 


+ +++ — 0.529820" x 1251)/0.510813 
= 1.84, with 5 degrees of freedom. 


; 
| 
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The individual ratios differ no more than would be expected by chance. 

After the manner of Section IX, formulae (13) and (14), an alter- 
native statistic for testing consistency is based on comparison of the 
ratios T/N. 


V(T/N) = + 4u + + + = 0.549534/N 
Hence 
X? = (194°/324 + 162°/268 + --- — 754°/1251)/0.549534 
1.33. 
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INTRODUCTION 


Methods for estimating variance components from unbalanced dat. 
are given in Henderson [1953] for both the random model and the 
mixed model. The calculations involved in the latter case are some- 
what tedious and usually not computationally feasible for data having 
many classes. This paper outlines the analysis for unbalanced data 
from a two-way classification with one classification considered fixed. 
Simplified computing procedures are presented suitable for a large 
number of levels in the random classification and a reasonable number 
of levels of the fixed classification. 


MODELS AND ESTIMATION 


The model for the two-way classification can be taken as 

= + By + + eine , 4) 
where 2z;;, is the observation and y is the general mean. aq; is the effect 
due to the 7’th level of the a-class, 8; the effect due to the j’th level of 
the 2-class, a@;; the interaction and e,;;, a random error term. We ‘vill 
suppose that the number of a-classes in the data is a, the number of 
B-classes is b, that there are n;; observations in the 7j’th subclass. and 
that s sub-clacses have observations in them. All terms except u ave 
taken as independent random variables in the random model with zero 
means and variances , 03 , and o2 respectively. The 6-ciassifica- 
tion wili be considered fixed in the mixed model, the variances invelved 
being ¢2 , 72, , and o? , the interaction terms being random for the 
particular set of fixed effects occurring in the data. In this case we are 
assuming an underlying multivariate distribution with 2, ;, having mean 


1Parts of this paper are based on the thesis submitted by the first author in partial fulfillment of 
the requirements of the Ph. D. degree. 
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value + 8; and variance + + : the covariance between 
Z,;, and for k k’ iso? + o%, and that between and 
forj jandk k’,iso?. 

Variance components can be estimated in both models by initia 
linear functions of sums of squares to their expected values, the sums of 
squares used being reductions in the total sum of squares due to fitting 
various elements of the model as if it were a fixed model. For example, 
fitting (u + a;) results in matrix equations of the form 


Pa = y, 
where y is the vector of a-class totals z;.. , and a is the vector of esti- 


mates of the a’s. P is non-singular, a diagonal matrix of order a, of 
terms n;.. Thus 


Pry, 
and the reduction in the ‘otal sum of squares due to fitting the a’s is 
R(u, a) = ay = y’P''y. 
This is the usual uncorrected sum of squares, namely 


Riu, a) = (2) 
Similarly 
Rw, B) = (3) 
R(u, a, B, aB) = , (4) 
= (5) 
and 


The variance components of the random model are estimated by 
equating the expected values to the observed values of differences among 
the above uncorrected sums of squares, namely (2)-—(5), (3)-(5), (4)- 
(2)-(3) + (5) and (6)-(4). These cannot be used in the mixed model 
because, apart from (6)-(4), their expectations then contain functions 
of the fixed effects. The within-subclasses sum of squares, (6)-(4), 
can still be used to estimate the error variance as in the random model, 
and ¢2 and o%,, can be estimated from the differences: 


Rix, a, B) — Riu, 


and 


Riu, a, B, a8) Rx, a, 8), 


| 
| 
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whose expectations are free of terms in the fixed effects. R(u, a, 8) is the 
reduction in the total sum of squares due to fitting a and 8 alone without 
the interaction terms. This reduction in sum of squares does not occur 
in the random model analysis and, when required for the mixed model 
analysis it usually invoives a considerable amount of computing. Its 
relationship to the random model analysis will now be shown and a 
simplified computing procedure derived. 


SIMPLIFIED EXPRESSION FOR R (yu, a, 8) 


Henderson ef al. [1959] have shown that fixed effects in a mixed 
model can be estimated by treating the random effects as fixed and 
maximizing the joint distribution function of the observations and the 
random effects. The resulting equations for the two-way classification 


without interaction are 
Q’ 


where P is a diagonal matrix, of order a, with terms n,. , R is a diagonal 
matrix of order b — 1 with termsn., , j = b, Q is a matrix of order a 
by b — 1 with terms n,; , j * b, y is a vector of the z,.. totals and z 
is a vector of the z_;. totals. R, Q and z’ have b — 1 columns because of 
omitting the equation for the last 8, appropriate to imposing the con- 
straint 8, = 0 in order to have a unique solution. R(u, a, 8) from these 
equations is 


a, 6) = 


= P (8) 


This expression for obtaining R(u, a, 8) requires inverting a matrix of 
order (a + 6 — 1), which is not feasible in many situations because a, 
the number of random classes in the data, is large. For example, the 
analysis of dairy production records in Searle and Henderson [1960] 
involved 688 herds (random) and 4 age groupings (fixed) so that a 
matrix of order 691 would need to be inverted for equation (8). How- 
ever, because of its special form, it can be reduced to inverting a 3 X 3 
matrix. Thus in the general case for the two-way classification the 
inversion of an (4 + b — 1) matrix can be reduced to inverting one of 
order (6 — 1) and this is a considerable reduction in the computing 
required especially when a, the number of random effects, is large, 
as is frequently the case. — 


y 
= (7 ) 
i 
7 
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The simplification of (8) proceeds as follows. Let 


D 
Then 

A=P'+ P'QDQ’P", 

B = —P'QD, (10) 


D = R 
pd I: 
=yP'y (11) 


The matrix P is easily inverted, being diagonal, and the only non- 
diagonal matrix to invert is D = (R — Q’P™’Q)™* of order b — 1. 

The first term in the above expression is R(u, a), and the second 
can be derived from equations (7). Eliminating 


a = — Q8), 
gives 

6 = Diz — Q’P’'y). (12) 
This is easily computed and can then be used to obtain R, , the re- 


duction in the sum of squares due to fitting the 6’s in this no-interaction 
model, as 


Ry = — (13) 
Thus from (11) R(u, a, 8) can be expressed as 
Riu, a, B) Ru, a) + Rs (14) 


the first term of which is part of the random model analysis and the 
second term comes from (12) and (13). 


EXPECTED VALUES 


The expectation of R(u, a, 8) can be found by using (8) and (9) 
and writing 


D 


= (y’Ay + z’By + y’Bz + z’Dz), 
= tr (Ayy’ + 2Bzy’ + Dzz’), 


: 
4 
4 
an 
| 
= | 
at 
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where tr(A) is the trace of the matrix A, the sum of its diagonal terms. 
After substituting from (1) for the vectors of totals y and z and also 
using (1) in (3), (4) and (6) it can be shown, with a certain amount of 
algebraic manipulation, that the differences used for estimating the 
variance components have expected values 


E(R(u, Qa, B) Riu, B)] 

= (n.. — + (kp — + (a — 
F{R(u, a, B, a8) Riu, a, 

= (n.. + (s —-a-b 
E{R(O) Ru, a, B, a8) ] 

(n.. — . 


The variances are estimated by equating the right-hand sides of 
these equations to the calculated values of the differences in the square 
brackets of the left-hand side. The term k, is one of the coefficients used 
in the random model analysis, namely k, = >>°., [}0%., n?,]/n.;.. The 
coefficient kg comes from the expected value of R(u, a, 8) and can be 
expressed in the form 


kg = 2tr (AU + 2BV + DW), (15) 
where U is a diagonal matrix of order a with terms ),°} n?, for 
i = 1--- a, V is a matrix of order a by b — 1 with terms ni, 
fort = 1 --- a,andj = 1--- b — 1, and W is a diagonal matrix of 


order b — 1 with terms n?, forj = 1. 


COMPUTING Rg AND kg 


R, is obtained from computing the following terms: 
(i) A square matrix C, of order (b — 1) with terms 


and 
= — j# = 1, eee , 


Computing c,; provides the check that }>°.., ¢;;- = 0, for all j, ne. 
row and column totals of the augmented C are zero. 
(ii) A column vector r, of order b — 1, whose terms are 


Dinid.., 


Computing r, provides the check that the sum of all the r,’s is ze: 


‘ 
. 
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(iii) A column vector, of order b — 1, of the estimates of the A's, 
6=C"'r. 
(iv) The “inner product” of the terms in the preceding two vectors is Rg: 
“b-1 
R, = 6'r = >> 6,7; = r’er. 
i=l 


This is equation (13), the C~* used here being identical to D. 

Expression (15) for ks is obtained using the elements of C~’, which 
we shall call d;;- , obtained in the calculating of Rg above. The special 
forms of the matrices involved in (15) enables kg to be written as 


where | 
A= 
= +. — 2ni5), 
and 


= Ai — — 
The f's, of which 3b(b — 1) in number are required, can be calculated 
for each 7 and summed, and the expression is well-suited to evaluation 
on a computer. Desk calculation can be arranged from the same 
formula when b is small, and the number of a-classes (random) is not 
too large. Calculating all of the $b(b + 1) /f’s enables the following 
check to be placed on them: 


b 
The f-values and }>; , can be obtained simultaneously with the terms 
of C and r; Rg is calculated as r’C~'r and kg is obtained using the elements 
of C~* as above. 


EXAMPLE 
The calculation of R, and kg is demonstrated for the small hypotheti- 
cal example shown in Table 1. 


The steps for obtaining R, are as follows: 
(i) The c-values are 


¢,, = 53 — 1°/10 — 37/20 — 127/40 — 127/50 — 25°/50 = 33.47, 
45 — 2°/10 — 6°/20 — 12°/50 — 25°/50 = 27.42, 


C22 


| 
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TABLE 1 
HyYpoTHETICAL EXAMPLE 
Number of observations n;; Mean 
observed 
j Totals values, 
1 1 2 3 4 nr x 
1 1 2 3 4 10 200 
2 3 6 4 7 20 300 
3 12 — 12 16 40 400 
4 12 12 13 13 50 500 
5 25 25 50 600 
Totals 
n.j 53 45 32 40 n.. = 170 
Totals of observed 
values, =. ;, 23000 23000 14000 19000 


C33 = 32 — 37/10 — 4°/20 — 127/40 — 137/50 = 23.32, 
= 40 — 47/10 — 77/20 — 16/40 — 13°/50 = 26.17, 
C2 = —1(2)/10 — 3(6)/20 — 12(12)/50 — 25(25)/50 = —16.48, 
C3 = —1(3)/10 — 3(4)/20 — 12(12)/40 — 12(13)/50 = — 7.62, 
and similarly 


C4 = —9.37, Ca = — 6.02, 
—10.78. 
The check on these values, namely >>°.., c;,, = 0, can be seen to hold 
true; for example 
Cir + Cio + C13 + Cry = 33.47 — 16.48 — 7.62 — 9.37 = 0. 


(ii) The r-values are obtained using the totals of the observations 
for the levels of the fixed effects and the means for the levels of the 
random effects: 


r, = 23000 — 1(200) — 3(300) — 12(400) — 12(500) — 25(600) = —3900, 


C23 = —4.92, C34 


r,= 23000 — 2(200) — 6(300) — 12(500) — 25(600) = —200, 
r,= 14000 — 3(200) — 4(300) — 12(400) — 13(500) = 900, 
r,= 19000 — 4(200) — 7(300) — 16(400) — 13(500) = 3200. 


A check on these is that their sum is zero. 
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(iii) 
33.47 —16.48 —7.62|" 
27.42 
~7.62 —4.92 23.32 
(iv) 
R, = r’Dr 


"053824 .036902 .025373 || —3900 | 
(—3900, — 200, 900) 036902 .063205 .025393 || —200 
| 025373 .025393 .056530 900 | 


736703. 


Calculating ks involves the elements of D and the f-values. The 
latter, together with the \,-terms are shown in Table 2. The calcula- 
tions for i = 1 are as follows. 


Ar = (1? + 2? + 3’)/10 = 1.40, 
fin = 17°14 + 10 — 2)/10 = .94, 
fioo = 2711.4 + 10 — 4)/10 = 2.96, 
hiss = 3°(1.4 + 10 — 6)/10 = 4.86, 
faz = 1(2)(1.4 — 1 — 2)/10 = —.32, 
fas = 1(3)(1.4 — 1 — 3)/10 = —.78. 


Similarly f,,.4 = —1.44, fi.23 = —2.16, fi.o4 = —3.68 and f,.34 = —6.72. 
The sums of these terms, over all values of 7, are shown in Table 2 
from which the f-values can be checked by the expression (16); for 
example 


505.8357 — 360.9718 — 113.1132 — 158.0607 

— 126.31 

—4°(1)/10 + 7°(3)/20 + 167(12)/40 + 137(12)/50 
‘= —(1.60 + 7.35 + 76.80 + 40.56). 

From the c’s and the f’s kg is now computed as 

ke = 45.79 + [505.8357(.053824) + 436.5532(.063205) 

+ 212.4332(.056530)] 
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TABLE 2 
Terms Usep IN CaLcuLaTING kg 


t fiu Siu 

1 1.40 .9400 2.9600 4.8600 5.4400 
2 3.05 7.6725 19.8900 12.0400 22.1725 
3 7.20 83 .5200 83.5200 7.2800 
4 9.14 101.2032 101.2082 112.0132 112.0132 
5 25.00 312.5000 312.5000 _ _ 


Total 45.79 505.8357 436.5532 212.4332 236.9057 


—fine —fias —fi.as —fi.ss 
1 .3200 .7800 1.4400 2.1600 3.6800 6.7200 
2 5.3550 2.3700 7.2975 8.3400 20.8950 11.1300 
3 60 . 4800 99 .8400 99 .8400 
4 42 .7968 49 .4832 49 .4832 49.4832 49.4832 56 . 9868 


Total 360.9718 113.1132 158.0607 59.9832 74.0582 174.6768 


— 2(360.9718(.036902) + 113.1132(.025373) 
+ 59.9832(.025393)] 
= 77.16. 


This procedure for obtaining Rg and kg may not appear greatly more 
straight-forward than inverting the matrix required in (8) which in 
this case is the 8 X 8 matrix 


20 3 6 4 

40 12 0 12 

50 12 12 13 

50 25 25 0} 

2 6 O 12 25 45 

4 12 13 O 32 


Advantages are apparent however, when one considers the case of a 


ay 
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large number of levels of the random classification 500 say, instead of 5. 
The matrix to be inverted would then be of order 503 while the pro- 
cedure outlined here would still only require inverting a 3 X 3. The 
calculation of its terms the c’s and of the terms for kg , the f’s, is still 
lengthy but can be accomplished separately for each 7 and summation 
made over i. This can be arranged quite straightforwardly for a desk 
calculator and is easily organized for an electronic computer such as 
the IBM 650, for example. 
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OPTIMUM SAMPLE SIZE IN ANIMAL 
DISEASE CONTROL’ 


A. W. Norpsxoe?, H. T. Davin H. B. 
Towa State University, Ames, Iowa, U. S. A. 


1. INTRODUCTION 


The objective of an animal disease control program would either 
be (a) to completely eradicate a disease or alternatively (b) to reduce 
the incidence of a disease to a low level and then to keep it in check. 
In the case of (b) the causative organism might not be completely 
stamped out. : 

If a disease has not obtained a strong foothold in a country, then 
objective (a) clearly is justified and perhaps the most drastic methods 
of control are in order. If, on the otherhand, a disease is widely dis- 
tributed both geographically and zoologically, then alternative (b) 
might be more realistic. This implies that the objective of a particular 
control program should take into account the cost of the program 
relative to the price of the consequence if there were no control. 

A case in point is the method of control for Pullorum disease 
(Salmonella pullorum) in poultry. The active form of the disease 
causes an enteritis in baby chicks which may result in high death 
losses. Chicks which recover may be carriers of the disease as adults. 
Female carriers may then transmit the organism via the egg to their 
progeny, where the disease again may become epidemic. Other species 
of domestic birds, and numerous wild species including pheasants, 
quail, and even foxes, cats, swine, cattle and man, have been reported 
as potential carriers of Pullorum. - 

Control of this disease is based on a blood testing technique. 
Samples of blood, collected from each member of an adult flock, are 
tested against an antigen of pullorum organisms for agglutinizing 
antibodies. Positive reactors to the blood test are judged to be disease 
carriers and therefore are removed from the flock and slaughtered. 

The blood agglutination test, carried out in most states for the past 
three or four decades, has now reduced the disease to one of relatively 
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minor importance. For example, the number of reactors reported by 
participants in the Iowa Control program was less than .07 of one 
percent in 1959 while in the Massachusetts program no positive reactors 
were reported in that year. 

When a disease such as Pullorum is widely distributed but at a low 
incidence on a state-wide or national level, it becomes problematical 
whether to continue past practices of completely blood testing all birds 
in a flock, or whether blood testing only a sample, or even no birds at 
all, of each flock might be more economical. 

The present paper investigates the extent to which 100 percent blood 
testing is justified, assuming random, that is binomial, infestation of 
adult flocks. Specifically, we attempt to determine the most economical 
flock proportion to be blood tested in any given year and locale, as a 
function of certain values and costs and of the prevailing binomial 
infestation rate. 

Economic optimizations in the presence of binomial a priori distri- 
butions have been studied previously, as for example in [1], [2], [3] 
and [4]. However, this prior work is almost entirely concerned with 
industrial inspection problems. The present paper is intended in part 
to illustrate the fact that such methodology, initially intended for 
industrial application, is equally applicable in the biological realm. 

We recognize that the assumption of binomial infestation might in 
certain instances be improved upon, by suitable contagion models for 
example. However, this assumption seems not unreasonable in the 
case of Pullorum; in this case, the adult flock is composed of birds 
which, though able to transmit the disease to their progeny through 
the egg, have low contagious influence as carriers. 


2. ALTERNATIVE FLOCK TESTING POLICIES; 
VALUES AND COSTS 


We consider the following one-parameter (the parameter is n) 
family of flock testing policies: 


Sample and blood test n birds out of a flock of N. If the 
sample contains no reactors (reaction assumed equiva- 
lent to infection), return the n birds to the flock and 
blood test no more. If the sample contains one or more 
reactors, blood test the entire flock; slaughter all dis- 
‘covered reactors. 


The relative economic worth of each of these N + 1 competing 
flock testing policies are now evaluated in the light of the following 
values and costs. 


i 
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i: blood test cost per bird, 

c: carcass value per bird, 

a(k): average value of a member of a flock of size N that 
contains k undiscovered reactors, 

b(k): value of a non-reactor from a completely tested 
flock of size N containing k discovered reactors. 


Note that a(0) = b(0) is the value of a member of a flock con- 
taining no reactors. Again, the value of a discovered reactor from a 
completely tested flock isc. Finally, the precise shapes of the function 
a(k) will depend on the likelihood that the progeny from a mating 
involving an undiscovered reactor will become acutely infected, thereby 
precipitating an epidemic in the progeny flock. The shape of the 
function b(k) will relate largely to loss of good-will. If good-will is 
not an important factor, it will not be unreasonable to set b(k) = 6(0) 
for all k; indeed, this is the assumption made at the end of Section 5. 


3. PRELIMINARY CONSIDERATIONS 
Computationally tractable forms of a(k) and b(k) are as follows. 


a(0) = A,a(k) =a for k #0, 
b(0) = A, b(k) = B for k #0. 


It is shown below that the most profitable n is either 0 or NV. If this 
can, for the moment, be assumed, then the derivation of the most 
profitable sample size (i.e. the choice between 0 and N) becomes a 
matter of simple arithmetic, at least for the limiting cases of very small 
or very large flock size N. 

Consider a large number M of birds, grouped into m flocks of 
N = M/mbirds. The case of very small flock size is typified by N = 1, 
while the case of very large flock size is typified by m = 1. 

For N = 1, the average value of the M birds (i.e. one-bird flocks) 
equals 


(1) 


M(1 — + Mra, if n=0, (2a) 
if n=N(=), (2b) 


so that 100 percent testing will be more (less) profitable than no testing 
according to whether 


— a) > (<) i. (3) 


For m = 1, the average value of the M birds (i.e. one flock of Af 
birds) equals 


: 
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M(1 — t)a + Mra = Ma, if n=0, (4a) 
M(1 — + Mae — Mi, if n= N(=M), (4b) 


so that 100 percent testing will be more (less) profitable than no testing 
according to whether 


B— a> (<) +1. (5) 


Expressions (2a) and (2b) arise from the fact that, under both zero 
and complete testing, a proportion (1 — 7) of the M one-bird flocks 
will be disease-free, contributing an amount M(1 — 7) A to average 
value. In addition, there will be a value contribution from the Mz 
one-bird flocks containing a reactor; this contribution will amount to 
Mra in the case of no testing and to Mzc in the case of 100 percent 
testing. Finally, there is the testing cost Mz that is incurred under 
100 percent testing. 

Expressions (4a) and (4b) arise from the fact that a very large 
flock will contain an approximate proportion 7 of reactors; hence, 
under 100 percent testing, there will be approximately Mz birds con- 
tributing carcass value c to the flock, and approximately M(1 — 7) 
birds contributing value 8. Similarly, if M is large enough to make 
Mrz large, the flock of M birds will contain at least one reactor with 
probability essentially 1, leading to an average per-bird value of a in 
the absence of testing. 

The criteria represented by (3) and (5) constitute an almost adequate 
solution for the case of the value assumptions given in (1). The further 
computations of Section 4 will serve only to validate the assertion that 
the most profitable sample size must either be zero or N, and will lead, 
as well, to the analogues of (3) and (5) for flocks of intermediate size. 


4. THE COMPUTATION OF THE MOST PROFITABLE SAMPLE SIZE 
FOR THE CASE OF CONSTANT VALUE DIFFERENCE 0(k) — a(k) 


This is the case typified by form (1) of the functions a(k) and b(k). 
The computations proceed as follows. Let 


¢ = the binomial probability of k reactors in a flock of 


= ) — (6) 


and let 


h,..» = the hypergeometric probability of obtaining n non- 
reactors when drawing a random sample of size n 
from a flock of N birds containing k reactors and 


i 
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(N — k) non-reactors = *) / ). (7) 


Then for any particular one of the (N + 1) alternative policies pre- 
sented in Section 2, say the policy corresponding to sample size n, 
Pr {number of reactors in the flock = k; number of re- 

actors in the sample = 0} = ¢,h,. , (8) 


and the “profit”? ensuing if the number of reactors in the flock equals 
k and the number of reactors in the sample equals 0 is 


— ni. (9) 


Hence the contribution to expected profit of the events involving no 
reactors in the sample equals the sum of the products of (8) and (9): 


N-n 


2 (Na(k) — ni] (10) 


k= 
Again, for the same sample size n, 
Pr {numbers of reactors in the flock = k; number of re-. 
actors in the sample > 0} = ¢-(1 — h,..), (11) 
and the “profit” ensuing if the number of reactors in the flock equals 
k and the number of reactors in the sample exceeds 0 is 
(N — k)-b(k) + ke — Ni. (12) 


Hence the contribution to expected profit of the events involving one 
or more reactors in the sample equals the sum of the products of (11) 
and (12): 


> ((N — k)-b(k) + ke — — (13) 


Total expected profit will equal the sum of (10) and (13), which, 
neglecting terms not involving n and using the fact that h,, = 0 for 
k > N — n+ 1, can be written 


N-n 


(N-[a(k) — + — ¢] + (N (14) 
Expression (14) is further reducible to 


(1 — m)"-{N-E[a(k) — b(k)] + E[k-(b(k) — c)] + (N — n)-4}, (15) 


where the expectation E[ ] is with respect to the chance variable k 
having a binomial distribution with parameters (V — n) and x. This 
type of expectation arises from the fact that 


Fra, 
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Further reduction leads, except for the additive term (1—1)*-N(8—a), 
to the expression 


(1 — x)"(T — nS) = P(r), (16) 

where 
T = Nia — B+ — + 4], (17) 
S=7(8—c)+i>0. (18) 


Differencing P(n) now yields 
A(n) = P(x + 1) — P(r) = (1 — — 1)S — x(T — S)], 


which shows that A(n) is negative for n < (T/S) + (l/r) — 1, and is 
positive for n > (7T/S) + (1/7) — 1. This implies that no P(n) for 
0 <n < N can be larger than the larger of P(O) and P(N), which in 
turn implies that the most “profitable” sample size n either is 0 or N. 
Establishing this fact (which is reminiscent of conclusions reached in 
[1] and [2]) was one of the two objectives set for this section in the last 
paragraph of Section 3. 

The second objective set for this section was to derive a condition 
analogous to (3) and (5) for determining the relative profitabilities of 
the two sample sizes 0 and N for intermediate flock size N. But this 
now is simply a matter of comparing P(0) and P(N), where P(n) is 
given by (16). This yields the conclusion that 100 percent testing is 
more (less) profitable than no testing at all according to whether 


(6 — a)[1 — (1 — > (<) eB (19) 


We note that (19) does indeed specialize to (3) and (5) for N equal, 
respectively, to 1 and to ~. 

Although details are outside the scope of this paper, it may be of 
interest to point out that condition (3) arises naturally in the computa- 
tion of the sequential Bayes test policy for the value assumptions (1). 
Consider the much (though no doubt impractically) enlarged set of 
flock testing policies consisting of all sequential plans with terminal 
acts A and R: 


A : Stop testing; collect ¢ for every reactor culled out so ‘ar; return 
non-reactors to the flock; eventually collect a or 8 per bird of 
the unculled portion of the flock, depending on whether or 
not this portion contains at least one reactor, unless the entire 
flock is reactor-free, in which case collect 1. 


: 
4 
| 


ANIMAL DISEASE CONTROL 623 


R : Test 100 percent; collect c for reactors, and 6 for non-reactors, 
unless the entire flock is reactor-free, in which case collect A. 

A straightforward application of the methodology given in [2] then 
shows that 100 percent testing is the most economic of the policies in 
this enlarged set if r(¢ — a) > i. However, the complementary pre- 
scription for no testing if r(¢c — a) < i does not apply in this case. 


5. THE COMPUTATION OF THE MOST PROFITABLE SAMPLE SIZE 
FOR THE CASE OF LINEARLY INCREASING VALUE 
DIFFERENCE b(k) — a(k). 


A value assumption alternative to (1) is 


ai) = OSESN, 
(20) 
bk) = A-~E(A- |, OSKEN. 


Using (15) [which was derived without reference to any specific form 
of a(k) and b(k)], expected profit now becomes, except for additive 
constants not involving n, 


pam = 1 +29), (21) 


where 
T = (A — — r)(N — 1) + — 0) + 71], (22) 
S=(A— >0. (23) 


The function P(n) given by (21) is best described in terms of the follow- 
ing five parametric cases. 

Case I: S > 0, T > 0. Asn increases from —o, P(n) rise steadily 
from — ©, crosses the n-axis at n = —T/S, equals T at n = 0, turns 
downward somewhere between n = —T7'/S and n = N, crosses the 
n-axis once more at n = N, turns upward somewhere beyond n = N, 
and approaches the n-axis asymptotically from below as n approaches 
+o. 

Case II: S > 0,0 > T > —NS. As n increases from — ©, P(n) 
rises steadily from —, equals 7 at n = O, crosses the n-axis 
at n = —T7/S, turns downward somewhere between n = —T'/S and 
n = N, crosses the n-axis once more at n = N, turns upward somewhere 
beyond n = N, and approaches the n-axis asymptotically from below 
as n approaches + ©. 

Case III: S > 0, T < —NS. As n increases from — ©, P(n) rises 
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steadily from — ~, equals T at n = 0, crosses the n-axis at n = N, 
turns downward somewhere between n = N and n = —T7’/S, crosses 
the n-axis once more at n = —7'/S, turns upward somewhere beyond 
n = —T/S, and approaches the n-axis asymptotically from below as n 
approaches + ©. 

Case 1V S = 0; T > 0. As n increases from —~, P(n) decreases 
steadily from + ©, equals 7 at n = 0, crosses the n-axis at n = N, 
turns upward somewhere beyond n = N, and approaches the n-axis 
asymptotically from below as n approaches + ~. 

Case V: S = 0, 7 < 0. As n increases from — ©, P(n) increases 
steadily from — ©, equals 7' at n = 0, crosses the n-axis at n = N, 
turns downward somewhere beyond n = N, and approaches the n-axis 
asymptotically from above as n approaches + ~. 

Cases IV and V are easily summarized as follows: If A = 6 (corre- 
sponding, as indicated at the end of Section 2, to the absence of the 
good-will factor), then the only policies in contention are 100 percent 
testing (n = N) and no testing (n = 0), and 100 percent testing will be 
more (less) profitable than no testing according to whether 


m(c — y) > (<)2. (24) 


It seems of interest to note the resemblance of (24) and (3). 

Cases III and V are summarized by: If T + NS < 0, it is most 
economical to test 100 percent. 

For Case II, the most economical sample size is the n between 
—T/S and N at which P(n) turns downward. This case is of interest 
since it establishes the possibility of an optimum sample size other 
than 0 or N. This possibility has already been noted in [3]. 

For Case I, the most economical sample size is either n = 0 or the 
n at which P(n) turns downward, depending on the relative magnitudes 
at these two sample sizes. Note that, for T = 0, the n at which P(n) 
turns downward must be greater than zero, so that, as in Case II, the 
most economical sample size will be other than 0 or N. 


6. EXAMPLE AND CONCLUSIONS. 


Consider the case when good-will is not a factor, i.e. the case of 
constant b(k). In this case both the formulation of Section 4 and that 
of Section 5 imply that only the two policies of zero and 100 percent 
testing are in contention, the choice between these depending on the 
direction of a simple inequality. 

Defining the critical testing cost i. to be the testing cost for which 
zero and 100 percent testing are equally economical, it seems of interest 
to compute 7, for both formulations, using comparable value figures. 
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A typical value for a non-reactor in a fully tested flock is $3, a 
typical average value of a member of an untested flock is $2.40, and 
typical values for c, r and N are $.50, 10°° and 500. 

‘Computing 7, in the spirit of Section 4, we therefore set A = 8 =$3 
and a = $2.40. Replacing 1 — (1 — x)” by Nx (allowable since Nx 
is small), criterion (19) then becomes $3 X 107° > (<) $2.5 X 107° + i, 
which means that 7, = $0.003, ie. that no testing is most profitable 
unless the cost of testing falls below 3 mills per bird. 

Computing 7, in the spirit of Section 5, we set A = 6 =$3. In 
addition, we interpret $2.40 to be the value of the linear function a(k) 
evaluated at k = Nz, the expected number of reactors in the flock. 
This implies a per-bird disaster value of y = —$6 X 10‘ for an un- 
suspected 100 percent infected flock. Criterion (24) then becomes 
$0.6 > (<)t, which means that 7. = $0.60, i.e. that 100 percent testing 
is most profitable unless the cost of testing rises above 60¢ per bird. 

In practice, the cost of testing is approximately seven cents per 
bird. Since this cost is well bracketed by the critica! costs 0.3¢ and 
60¢ derived above, we learn that the shape of the value function a(k) 
fand of course also that of b(k)] must be determined rather accurately 
if the methodology presented here is to be applied. 
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NUMERICAL ASPECTS OF THE REGRESSION 
OF OFFSPRING ON PARENT" 


H. E. McKeEan anp B. B. BouREN 


Population Genetics Institute, Purdue University 
Lafayette, Indiana, U.S. A. 


INTRODUCTION 

In an earlier paper (Bohren, McKean, and Yamada, [1961]) three 
currently employed techniques for estimating the regression of off- 
spring on parent, and thereby heritability in the narrow sense, were 
compared and contrasted with respect to their efficiencies of estimation. 
The general conclusion, based on theoretical considerations and an 
empirical study of five generations of a closed poultry flock (Yamada, 
Bohren, and Crittenden, [1957]), was that under the circumstances 
considered, the method of regression of offspring means on parent’s 
record (method 1) was inferior to the method of regression of individual 
offspring on parent (method 2) (with the parent’s record repeated once 
for each of its offspring) and to (method 3) the Kempthorne-Tandon 
technique (Kempthorne and Tandon, [1953)). 

The success of the Kempthorne-Tandon technique depends upon 
knowledge of a parameter p, the correlation between deviations of two 
offspring of the same parent from the predicted breeding value of the 
parent, and its expected superiority over the second method depends 
upon the magnitude of p. Usually p is guessed in the light of prior 
knowledge, and weights are assigned to the families according to the 
guessed value of p. In the first paper it was shown that, under the 
assumption of all genetic variance being additive, p < .067 or < .079, 
depending upon whether the mating structure is random or hierarchal. 

The results obtained in the previous paper specifically depended 
upon the particular distribution of family sizes encountered in the five 
analyses, and upon the accuracy of the estimated values of T=p/(1—p). 
The purposes of this paper are to consider the efficiency loss incurred 
by mis-guessing p in the Kempthorne-Tandon technique, and to in- 


vestigate the factors involved in the relative efficiency of the other 
two methods. 


THEORY 


We consider a breeding experiment in which s sires are selected 
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from the population, sire 7 being mated to a random sample of d, distinct 
dams. The mating structure is then hierarchal, with the progeny of 
sire i having phenotypic values given by the model 


Yin =u + (1) 


which is equation (11) of the previous paper. 

It has been pointed out that the three methods under consideration 
are merely special cases of the general unbiased weighted regression 
estimation procedure. The difference between the methods involves 
only the weights applied to each progeny-group deviation: (1) w,; = 1 
for the progeny means on parents technique, (2) w,; = n,; for the 
repeated parents technique, and (3) w,; = ,;/(1 + m4;7) for the 
Kempthorne-Tandon technique, where 7 is a guessed value of T. If 
7 = T,, the third technique is the minimum variance technique, whereas, 
if r = 0, the third technique reduces to the second. Furthermore, 
when the family sizes are equal (n;, = nig = --- , for all 2), all three 
methods are identical in efficiency. Since most experimental data 
will involve unequal family sizes, interest will center on considering 
this situation. 

The question, “When may I use method 1 with little loss in effici- 
ency?”, is a pertinent one. First, let us examine the optimum choice 
of weights for which the variance of the estimate of 8 will be minimized. 
These optimal weights w* (say), where w* = n,;/(1 + n,;T), will be 
approximately equal (hence method 1 appropriate), for JT’ > 0, under 
one or more of three distinct sets of circumstances: 

1) All n,; are large. This follows immediately from 


tan 

LL + T 
2) S? , the variance of the family sizes, is zero or very small. 
3) T is large. This follows from the fact that 


* 
(or p>1) 
independent of n,;; and n,-;- ; thus for large T, w* = w*;.. In view of 


this, method 1 may also be considered as a special case of method 3 
where 7, the guessed value of 7’, is allowed to approach ~. 

It is of interest, therefore, to determine a readily accessible statistic 
which depends upon 7 (the average dam family size), S? , and 7’, upon 
which a decision to use or not use method 1 as opposed to method 2 
may be based. 

It is easy to show that the coefficient of variation among the optimal 
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weights w is approximately 


S, _ C.VAn) 
where C.V.(n) = S,/ni = coefficient of variation of the dam family 
sizes. Since ¥(T) is a function of all three of the statistics mentioned 
above, it is of interest to consider the relationship between relative 
efficiency (Min V(8)/V(§)) of the various techniques and ¥(T) [even 
though ¥(T) may seriously differ from the true coefficient of variation 
of the 
For any choice of r, we may define 


= C.V(n)/(1 + fiz). (2) 


Clearly, ¥(0) = C.V.(n), ¥(©) = 0, and ¥(r) is a monotonically de- 
creasing function of 7. Since the range 0 < +r < © corresponds one- 
to-one to the range in y of C.V.(n) > yw > 0, it is more convenient to 
express graphically a relationship between y, 7’, and relative efficiency 
rather than 7, 7’, and relative efficiency. 

Given a fixed set of family sizes n,, , M2, --- the relative efficiency 
is a function of 7, the true unknown parameter of the population, and 
t, the guessed value. Because of the one-to-one relationship between r 
and 7(7), the relative efficiency may be expressed as a function of 7 
and y. The explicit form of this function may be obtained by solving 
equation (2) for 7 and substituting the resulting expression in y into the 
equation for relative efficiency: 


where £; = estimated regression within the ith sire using r as the guessed 
value of 7’, 

6* = estimated regression within the 7th sire using r = 7' (optimum), 
and 


1+n.,T — x.) 


ot = — p)/D wi(X,, — 


which correspond to equations (15) and (16) of the previous paper. 
The graph of a typical relationship between y, 7’, and relative 

efficiency is given in Figure 1. The z- and y-coordinates of any point 

on the surface pictured represent respectively, the coefficient of varia- 
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FIGURE 1. 


GRAPHICAL REPRESENTATION OF A TYPICAL RELATIONSHIP BETWEEN z = y¥, y = T, 


AND z = RELATIVE EFFICIENCY. 


tion of the weights resulting from a choice of 7, and the true value of 
the parameter 7, whereas the z-coordinate is the resultant relative 
efficiency. Several interesting properties should be noted: 


1) 


2) 
3) 


+ 


— 


5) 


6) 


The points in the (x, y)-plane corresponding to s = T (optimal 
efficiency) lie on a curve whose equation is x = ¥(T). The 
surface representing relative efficiency thus has a ridge which 
coincides with the projection of this curve onto the plane z = 1. 
The relative efficiency falls off from this ridge in all directions. 
The relative efficiency of methoa 1(r = ©) is a strictly in- 
creasing function of T. 

The relative efficiency of method 2(r = 0) is a strictly decreasing 
function of T. 

The relative efficiency of method 3(0 < +r < ©) is strictly in- 
creasing for 0 < < 7 and strictly decreasing for r < T < 
For any fixed value 7’, relative efficiency is a strictly increasing 
function of 7 for 0 < +r < T and strictly decreasing for r > T. 
This family of curves, generated by all possible 7 values, ex- 
presses graphically the consequences of mis-guessing 7’. 
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7) If C.V.(m) is small, the whole’ picture is condensed along the 
z-axis so that mis-guessing the value of 7 will not result in 
serious loss in efficiency. If C.V.(n) = 0, the surface is degenerate, 
and reduces to the line (x = 0, z = 1). Alternatively large 
values of C.V.(n) may result in serious efficiency loss if | 7 — T | 


is large. 
AN ILLUSTRATION 


In order to consider the effect of a poor guess for p, we will consider 
the poultry population previously studied by Bohren, McKean and 
Yamada [1961]. Five years data, 1952-1956 inclusive, are included 
in the analysis. The information pertinent to the population size and 
structure are given in Table 1 along with the average family sizes and 
the coefficients of variation of n. The relative efficiencies of the three 
methods for the year 1952 are shown in Table 2 for a range of values 


TABLE 1 

Tue Famity STRUCTURE OF THE POPULATIONS STUDIED. 

No. No. No. Average C.V. (Family 
Year Sires Dams Pullets Family Size Size) 
1952 10 106 605 5. .74 
1953 10 92 877 9.5 43 
1954 1l 78 573 7:33 Al 
1955 18 108 852 7.9 47 
1956 20 132 715 5.4 .54 


of T and 7 between zero and .07. This covers most of the theoretical 
range of p, under the assumption of additive genetic variance, as found 
by Bohren, McKean, and Yamada [1961] and includes the estimates of 
T obtained from these data. It is seen that in no case is the efficiency 
of either the repeated parent technique (r = 0) or the Kempthorne- 
Tandon technique less than .97. If an intermediate value of 7 is chosen 
under the Kempthorne-Tandon technique (say .03) it is seen that the 
efficiency does not decline below .99. The efficiency of the means 
technique is seen to be consistently poorer than the other two techniques. 

Analogous tables were generated for the years 1953-1956. Es- 
sentially identical results were obtained for all five years analyzed. 
In all years the relative efficiency of the Kempthorne-Tandon technique 
remains above 97 percent for T < .07 andr < .07. 

Method 1 was consistently less efficient than method 2 in all years 
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TABLE 2 
RELATIVE EFFICIENCIES OF THE THREE METHODS OF EsTIMATING 8 
RELATIVE TO THE MINIMUM VARIANCE ESTIMATOR OF A(T’ = 1) 
FoR VARIOUS SMALL VALUES OF T' FoR 1952 Data. 


Method 1 Values of 7 used in the KT technique (Method 3) 


Tech.) .00? 01 02 03 - .04 05 06 07 


27 = 0 is Method 2 


over the range 0 < T < .07, but was relatively more efficient than. in 
1952. Letting EF denote the efficiency of method 1 (means), then, as T 
rises from zero to .07, in: 


1952, E rises from .5836 to .7126, 
1953, E rises from .9142 to .9615, 
1954, EF rises from .9046 to .9587, 
1955, E rises from .9370 to .9682, and 
1956, E rises from .8058 to .8703. 


Table 3 gives a direct comparison of the relative efficiencies of 
methods 1 and 2 when large values of T are involved, which may, but 
do not necessarily, arise when maternal effects, non-additive genetic 
effects, and/or environmental correlations are present. It is interesting 
to note that the T-value at which the two methods are equally efficient 
varies considerably due to differences in the distribution of dam family 
sizes, varying from approximately T = .09 in 1953 to T = .32 in 1952. 
The corresponding values of ¥(7), the coefficients of variation of the 
optimum weights w* , show much more year-to-year stability, how- 
ever, ranging from a low of .23 in 1953 to a high of .27 in 1952. 

Let 7, represent the value of T for which the two methods are 
equally efficient. The values ¥(7,) will not ordinarily be limited to as 
small a range as observed here (see Table 4). It is quite possible to 
have much lower values of C.V.(n), which is an upper bound for ¥(T7%). 
Hence one cannot claim any invariance properties for ¥(7'>). On the 


(Means 

00 | .5836 | 1.0000 .9988 .9958 .9919 .9872 .9822 .9770 .9716 

01 | .6083 .9988 1.0000 .9991 .9969 .9938 .9901 .9861 .9818 

02 |  .6302 .9991 1.0000 .9993 .9976  .9951 .9922 .9888 

03 | .6499 .9969 .9993 1.0000 .9995 .9981 .9961 .9936 

04| .6678 .9937 .9976 .9995 1.0000 .9996 .9985 .9968 4 

05 |  .6840 9819 .9901 .9951 .9981 .9996 1.0000 .9997 .9987 

06 | .6989 9766 .9860 .9922 .9961  .9984  .9996 1.0000 .9997 

07 .7126 9713 9818 .9889 .9937 .9968  .9987 .9998 1.0000 

: 
S. 
d. 
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TABLE 3 
RELATIVE EFFICIENCIES OF METHODS 1 AND 2 FOR LARGE VALUES OF T’. 
Year Method .10 -20 .30 .40 .50 
1952 1(Means) .748 .825 .870 .900 .920 
= 0) .956 .912 .880 .856 .838 
1953 1( Means) .970 .984 .990 .993 .995 
2r = 0) .963 .932 .912 .898 .889 
1954 1(Means) .968 .984 .990 .993 .995 
2r = 0) .971 .947 .932 .921 .914 
1955 1(Means) .974 .985 .990 .993 .995 
2(r = 0) .974 .952 .937 .926 .918 
1956 1( Means) .887 .924 .945 .958 .966 
2(7 = 0) .975 .945 .920 .901 .885 
TABLE 4 
ILLUSTRATION OF THE UTILiTy oF ¥(7') IN CHOOSING 
BETWEEN METHODS 1 AND 2. 
Estimated Standard Error 
Year f KT) (To) Method 1 Method 2 
1952 059 55 .121 -102 
1953 029 34 .23 .085 .083 
1954 051 .30 .24 .152 .148 
1955 055 33 -26 .086 .084 
1956 034 46 .24 .108 .094 


basis of information available the recommended procedure for selecting 


an estimation technique is as follows: 


A. From previous experience, or by any other reasoning, assume a 
value 7 for T. 
B. If possible, use the Kempthorne-Tandon technique, constructing 
weights using this assumed value r. 

C. If it is neither possible nor convenient to use method 3, compute 
¥(r). For large values of ¥(r) (say > .4) preference is for method 
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2 whereas for small values of ¥(r) (say < .1) method 1 is preferred. 
For intermediate values of ¥(r) there should be no substantial 
differences between methods one and two. 


To illustrate the use of this approach, consider Table 4. Here, 
estimates of 7’, obtained by Bohren et al. [1961], for each of the years 
1952-1956, are used to calculate ¥(T). For comparison, ¥(7') is also 
entered, together with the estimated standard errors of 8. The values 
¥(T.) indicate a preference for method 2 in each case, but 7, is obtained 
only after tedious calculation. The rule of thumb suggested would 
lead to strong preference for method 2 in 1952 and 1956 and no preference 
in the other years. 

The results of these investigations indicate that in many circum- 
stances a poor choice of technique can result in serious loss in efficiency 
in the estimation of 8, and thus of heritability in the narrow sense 
(h? = 28). 

It would seem that, once the expense of experimentation and data 
collecticn has been incurred, the investigator is under an obligation to 
obtain highly efficient estimates. However, in some instances the loss 
in efficiency is negligible if one of the two standard techniques is used. 
The use of ¥(r) to aid in the choice between the standard methods 
seems promising. 
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SOME HYPOTHESES CONCERNING TWO 
PHASE REGRESSION LINES 


P. SPRENT 


East Malling Research Station, 
Maidstone, Kent, England. 


1. INTRODUCTION 


The regression of a growth measurement y on time z can often be 
reasonably represented by two intersecting straight lines, one being 
appropriate when z takes values below and the other when z takes 
values above a certain fixed but often unknown value corresponding 
to the intersection. Such regressions are here called two-phase re- 
gressions, the intersection of the phases being referred to as the change- 
over point, and the value of zx at which it occurs being called the change- 
over value. 

Situations in which such regressions might occur include the onset 
of a disease resulting in a reduced growth rate; the application of a 
treatment having an immediate stimulating or inhibiting effect; the 
occurrence of an extremely hot or cold day or some other change in 
external conditions; physical injury of an organism. 

In a study of the compatibility of peach scions on plum rootstocks 
Garner and Hammond [1938] noted that the peach variety Hale’s 
Early developed at constant but different rates on compatible and 
incompatible rootstock-scion unions up to a certain date. After that 
date the growth rate in the compatible case continued at a new constant 
rate, whilst in the incompatible case all growth then ceased. 

Thus for a compatible union there was a typical two-phase regression, 
whilst for the incompatible union a rather special case occurred in 
which the slope of the second phase was zero. A further example is 
given in Section 4 in which the date of phase change in relation to 
time elapsed after application of treatments is of interest. 

If x and y are-growth measurements on two different parts of the 
same organism, and Huxley’s allometric growth law operates (i.e. there 
is a linear relation between log z and log y) it has sometimes been 
found that sudden changes in slope occur. Reeve and Huxley [1945] 
have discussed this situation with reference to changes in growth 
equilibrium in crustacea at sexual maturity. Skellam et al. [1959] also 
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refer to such phase changes in allometric growth, mentioning earlier 
work of Teissier on this subject. 

A biologist will often postulate a two-phase linear regression rather 
than some alternative such as a parabolic one on largely intuitive 
grounds, and his decision on this point must be to some extent a matter 
of experience and common sense, as is generally the case in selecting 
appropriate hypotheses and models for statistical examination. In 
many cases a two-phase regression can only be a reasonable approxi- 
mation, adequate for many purposes, but by no means a complete 
description of what is taking place. 

It is likely that two-phase regressions also arise in economic and in 
industrial production problems. 


2. HYPOTHESES 


In the present context questions of the following type are likely to 
interest the experimenter. 

A. Necessity for, and adequacy of, a two-phase fit. 

B. Investigation of change-over points and values for individual 
regressions. 

C. Relationships between lines within each phase for several re- 
gressions. 

D. Relative positions of change-over points for several regressions. 
For example, interest may attach to one or more of the following 
hypotheses. 

In relation to A: 

AI. One phase adequate (i.e. a conventional linear regression 
suffices), 

AIT. A two-phase linear regression is adequate; 

In relation te B: ; 

BI. The observed change-over value is consistent with a specified 

theoretical value x = 

In relation to C: 
CI. The lines within a phase are all parallel to one another, 
CII. The lines within a phase all coincide; 

In relation to D: 

DI. The change-over points are collinear, all lying on a specified 
line y = a + (or = if infinite), 

DIT. As for DI except that a (or y) not specified, 

DIII. The change-over points are collinear on an unspecified line. 

In practice, anyone accepting AJ would proceed no further in this 
contest. If AZ were rejected and AJJ accepted, a test of BI might be 
required, and with more than one set of data an experimenter may wish 
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to proceed to some of the C or D hypotheses. For instance, he might 
test and accept CJ, reject CII as the result of a test (or without test if he 
considered it incompatible with his data) and then test and accept 
DIT for the case of infinite slope. If x were a time measure, acceptance of 
DII with 8 infinite weuld imply that all phase changes occured at the 
same time, and is thus likely to be of practical interest. If one wants 
to know if a phase change can be associated with an event occurring 
at a particular x value, such as the application of a treatment, or the 
occurrence of a very hot day, a test of BI is required. 

The above hypotheses are by no means exhaustive. Concurrence 
within a phase, for instance, may be of more interest than parallelism. 

In this paper two assumptions that will be made are (i) the classical 
regression assumption that within each phase the observed y are normally 
distributed with mean zero and standard deviation o about the true 
regression line, and (ii) that it is known between which two observed x 
the change-over occurs, but not precisely where between these two 
values it occurs. 

This latter assumption calls for some comment. As will be seen, 
it leads to certain simple tests, and in many practical cases it is not 
unreasonable. Sometimes observations can only be taken at widely 
spaced 2x-values, and it may be clear that the change-over has occurred 
between two such observations, without a more precise statement 
being immediately possible. 

A difficulty arises if lines are fitted to each phase making the as- 
sumption (ii), and the intersection of the fitted lines has an x co-ordinate 
falling outside the assumed range. The experimenter must then decide 
whether to attribute this to sampling errors, or an incorrect assumption, 
and determine the appropriate procedure accordingly. This matter is 
not pursued further here. 

If the second assumption cannot be made because of doubt about 
the allocation of a small proportion of points near the change-over, 
the author conjectures that the methods of Quandt [1958] could be used 
to allocate these doubtful points without seriously upsetting the pro- 
cedures suggested here. Quandt deals with the allocation problem for 
a two-regime y-on-z regression where the x occur in a known order 
with respect to a time measure ¢, and the phase change occurs when a 
certain t-value is reached. In the present case x = t. 

This paper deals mainly with questions classified under B and D, 
but we first briefly discuss hypotheses involving A and C. 

If one has an independent estimate of error variance, the test of A/ 
is simple and standard, being a test for significance of deviations from 
regression. If deviations from regression provide the only estimate of 
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error variance and indicate significance of linear regression, some judge- 
ment is needed to decide whether a two-phase regression (or any other 
fit such as a quadratic regression for example) is worth considering. 
The test of AJ is likewise essentially a test of the adequacy of two 
linear regressions, one in each phase. 

If AJI is accepted for more than one set of data and tests for 
parallelism or coincidence (C7, CII) are required, these may be per- 
formed by standard methods for each phase separately. 

We now consider other tests. 


3. TESTS 


Throughout summations over r run from 1 to 2 (corresponding 
to the two phases); those involving s run from 1 to n (corresponding 
to n sets of two-phase regressions); those involving k run from 1 to 
m,, for the appropriate set and phase, m,, being the number of observa- 
tions in the rth phase of the sth set. Where it is obvious how many 
equations (usually 2, n or 2n) of a given form exist, only typical ones 
are given. 

The first test considered is that of BJ, that is compatibility with a 
fixed change-over value, assuming AJJ. The test is based on the 
comparison of the fit by standard regression techniques of uncon- 
strained lines in each phase, with that obtained if the lines are con- 
strained to intersect at a point whose x co-ordinate y is given. 

Appropriate models for the first and second phases respectively are 
E(y) = a, + 8,2 and E(y) = a, + 6.2, where the constraint to give an 
intersection at = y — a, + — B,) = 0. 

In view of the regression assumptions the method of least squares 
is appropriate, and a, , a. , b, , b2 , the estimators of a, , a2 , 8; , Be, 
are to be determined so as to minimize 


subject to 


a, — a, + y(b, — 6,) = 0. (1) 


Although not essential, introduction of a Lagrange multiplier » 
preserves a certain symmetry, and gives rise to many expressions 
similar to familiar regression results, together with adjustment terms. 
With such a multiplier the function to be minimized becomes 


a, — btn)? + — a, + — 


whence, computing the partial derivatives, normal equations are ob- 
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tained which can easily be reduced to the form 
DX Yate — 4,2. — b,27,) + dy = 0, (2,3) 


a, = y,. — b.z,, + (-1)""A/m, , (4,5) 


where y,. = > Y-x/m, , etc., m, being the number of observations in 
the rth phase. The remaining normal equation is (1). Subtracting (4) 
from (5) and using (1) gives 


= wld + y) — ¥)], 


where w = mm./(m, + m2) = mm./m, say, andd = y.. — . 

Substituting for a, , a2 , \ as given by (4) to (6) in (2) and (8) gives, 
upon writing Cz,y, , Cx? , etc., for the ordinary corrected sums of 
products and squares, 


(e.g. Cxi = (tu — ete.): 


+ — — — — y)be 
= Cry: w = ¥); 


— w(z,. — y)(t2. — + + — 7)" 
= Cryo + — 7). 


Having solved these equations for b, and b, it is easy to obtain 
values for a, , a2 , \ if required from (4), (5) and (6). This constrained 
fit is then compared with the unconstrained fit for which of course, 
the slopes b’ , (r = 1, 2) are given by b’ = Cz,y,/Cz?. 

If there are m = m, + mz, observations, the sum of squares for 
departures from unconstrained regressions has m — 4 degrees of freedom, 
two regression coefficients and two means having been estimated. For 
the constrained fit, the corresponding sum of squares has m — 3 degrees 
of freedom since only three independent parameters have been esti- 
mated in view of (1). It is easily shown that the difference, with one 
degrees of freedom, has a sum of squares given by 


Dd — b,)Cz,y, + dd. (7) 


This sum of squares may also be obtained by subtracting 
2 b,Cz,y, — — 


from the unconstrained regression sum of squares. It is worth noting 
that the constrained lines will no longer in general pass through their re- 
spective centres of mass. 
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An F-test of (7) against an appropriate error variance is required. 
A significant value indicates a worthwhile improvement in fit if the 
constraint is removed, and in that event BZ would normally be re- 
jected for the 7 under consideration. 

The above test may be extended if there are n sets of data, and 
given change-over values y, , --- , Y, postulated, each one associated 
with a given set. The regression coefficients may be estimated separately 
for each set by the procedure above, using the appropriate constraint 
each time. If desired, a joint test of the fit compared to unconstrained 
fits for all sets within each phase may be made. The sum of squares 
(7) is computed for each set, and these are added together to give a sum 
of squares with n degrees of freedom which can be tested against an 
appropriate error variance. If y, = --- = y, it will be noted that this 
test is essentially a test of DI with 8 infinite. 

Turning now to DIJI there is little difficulty in developing a test 
if the change-overs are assumed to lie on a line of infinite slope (i.e. 
to occur at the same unknown z-value) providing AZZ and CI are 
accepted (i.e. all lines within each phase are parallel to one another). 
It will be assumed that all lines in the first phase have the same slope 
b, , and all lines in the second phase have the same slope by , (b,¥ b;). 
A comparison has to be made between a fit constrained so that all 
change-over points lie on a line x = y’, and a fit not so constrained. 
(Note y’ is not specified, and must be estimated.) 

The model now is 


=a, + Bx and Ely) = az, + Box 


for the respective phases, where s = 1, --- , n, and there are n con- 
straints which may be written 
Q2, — a, + — = 0. (8) 


As 8, , 82 and y’ are independent of s, it is convenient to replace (8) 
by the n constraints 


Qo, Ais = 


the geometrical significance of which, in relation to intercepts, is obvious. 
The function to be minimized in order to estimate these parameters 
(including ¥ in this case) is 


L, = > — Gr, — +2 — — 


the A, being Lagrange multipliers. Differentiation gives rise to a set 
of normal equations, e.g. 
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ars = Yre. + > 


>. = 0, 


— a, —c=0. (9) 


These lead, after some straight-forward but tedious algebra to the 
equations 


[do + — = — Com. d, (11) 
+ Cz. + > + CoB. d, (12) 


where Cz’, , Cz,,y,, are the usual corrected sums of squares and pro- 
ducts for the sth line of the rth phase, and 


W, = + Me,), 

Cote = — ote. (Do 
Cun. d= d. — w, d,), 

d, = You. — Yre., Cte. 


In the above all summations are over s. 

It is easiest to solve (11) and (12) for b, and b, , then obtain c from 
(10). For testing purposes it is convenient to know the \, , which can be 
obtained from the easily deducible relations 


A, = — bet... + — (13) 


Putting m = >>, >>, m,, , the total number of observations, then 
for unconstrained parallel lines the sum of squares for deviations from 
regression has m — 2n — 2 degrees of freedom. When the constraints 
(8) apply we fit only n + 3 independent parameters and thus have 
m — n — 83 degrees of freedom for deviations from regression. The 
difference between these two deviation sums of squares has n — 1 
degrees of freedom, and is obtained as 


Dd (bt — + rAd, , (14) 
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where b/ , bi are the regression coefficients for unconstrained parallel 
lines. 


Equations (11) and (12) simplify somewhat if the z’s are the same 
for all sets of data. 


4. EXAMPLE 


R. M. Fulford (personal communication) has made counts of the 
number of nodes in terminal buds of spurs of the apple variety Miller's 
seedling from the first or outermost budscale to the youngest pri- 
mordium at intervals after the application of three defoliation treat- 
ments. The treatments were applied on different dates, and the counts 
were made at intervals of seven or more days throughout the 1957 
growing season. On each counting date at least six buds from trees 
that had received a given treatment were examined. The counting 
dates were not always the same for all treatments, but were generally 
close together and covered a comparable time span. 

The regression of y, the number of nodes, upon z, the number of 
days after May 24th, was found to be consistent with the hypothesis 
of three sets of two-phase parallel regressions (ie. AJJ and Ci were 
satisfied). This was tested by standard methods which are not re- 
produced here. It appeared from the data that the change-overs might 
well have occurred on the same calendar date, rather than at, say, the 
same fixed number of days after the application of each treatment. 
In each phase of each regression there were observations for between 
six and eleven distinct counting dates x. 

The raw data are not reproduced here; but relevant means, sums 
of squares and products, weights, etc., are given. 

The content of Table 1 was computed directly from the data, and 
that of Table 2 from that of Table 1. 


The coefficients in equations (10) to (12) are easily got from these 
tables, and give 


c = 8-260 — 97-856b. + 39-578), , 
56926-21b, — 777-54b, = 11388-08, 
—777-54b, + 19057-96b. = 954-93, 
whence b, = 0.2008, b. = 0.0583, c = 10.502. From (13) we then get 
A, = —4-404, Ae = —5-61l, As = 10-017. 


The d, , apart from rounding errors, sum to zero, providing a useful 
check. Computing (14) we find it comes to 7.32. This sum of squares 
has 2 degrees of freedom. An appropriate error mean square is provided 
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by that between counts on the same date within treatments and has 
the value 1.08 with 230 degrees of freedom. This routine computation 
is not shown. The appropriate F is thus 3.66/1.08 = 3.39, and this is 
significant at the 0.05 probability level. Thus there is evidence for 
rejecting the hypothesis of a common change-over date, although it 
was clear from the data that the change-overs occurred nearer to a 
common calendar date than they did to a constant time lapse after 
application of each treatment. In the computations, particularly in 
calculating (14), care is needed to avoid excessive rounding errors. 


5. RELATION TO FIELLER’S THEOREM 


It is evident that a test of consistency of the observed change-over 
value with a specified x = y, (hypothesis BJ) is equivalent to a test 
about the intersection of two regression lines, and the test developed 
in Section 3 might be expected to be related to Fieller’s [1940, 1944] 
theorem on fiducial limits of a ratio. Fieller has shown that the fiducial 
limits for 6 = y/zx, where the errors in y and z are normally distributed 
with zero mean, are given by the roots of the quadratic in 6 


(y? — — 20(zy — + — #,.) = 0, 


where v,, , Vz , Uz, are estimates of the variances and covariance of y 
and z, each based on n degrees of freedom, and ¢ is the appropriate 
value of Student’s ¢ for the required fiducial probability with n degrees 
of freedom. 

Applying the theorem in the situation envisaged in BI, Fieller’s 
y = a, — af, andhisz = bf — bj, in terms of the parameters of un- 
constrained regression lines. The estimated variances and covariances 
of these, using an appropriate error mean square s° say, are easily 
obtained. On the criterion of Section 3 for testing BJ, a chosen y will 
just be significant if 


(bf — b,)Cz,y, + rAd — #s? = 0. (15) 
From Fieller’s theorem it follows that such a ¥ is given as a root of 
the quadratic in y, 
[y(b{ — bs) — (ah — — 


[E+ Cat (16) 


The left-hand side of (15) can be shown to be a quadratic in +, 
and thus (15) and (16) will be equivalent if it can be shown that their 
left-hand sides are identical for three distinct values of y. This is 
easily done when y = 2, or 22. or (ai — a{)/(b{ — bj). 
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Thus using Fieller’s theorem to establish fiducial limits for y pro- 
vides a set of y's (i.e. those between the limits) which would be accep- 
table under the hypothesis test put forward here. 

With either Fieller’s theorem, or the approach used here, difficulties 
arise if y is infinite, but in practice this causes no trouble in the present 
context as an infinite ~ would imply the identity of the two phases, 
i.e. acceptance of AJ, or alternatively a pair of parallel regressions, 
rather than the abrupt change of slope implied by 7 finite. 

Whilst Fieller’s theorem does not link directly with the test already 
considered for a common unspecified change-over date for all of n sets 
of data, it does provide a solution to a closely related problem, namely: 

Assuming there is a common-change over date, what are the p 
percent fiducial limits for this? Equations (10) to (12) provide 
estimates c, b, , bs , of vy, 8: , 8. respectively. The 7’ of (8) is related 
to y by the expression 


(8, 
thus y’ may be estimated as 
ce’ = c/(b, — be). (17) 


If the value of c given by equation (10) is put in (17) and the variances 
and covariance of numerator and denominator obtained in the usual 
way, Fieller’s theorem may be applied. 


6. DISCUSSION 


An important characteristic of the tests derived in Section 3 is that 
they involve linear constraints upon the parameters. This would not 
be the case for instance if one wished to test DJJ with £8 infinite but 
without the restriction to parallelism introduced in Section 3. Then 
the n constraints would become 


G2, — a,, = c(b,, — bz,), 


all terms being unknown. Further, if one wished to test DJJ with B 


finite, and no restriction to parallelism, the constraints would take the 
form 


These do not easily reduce to a workable system of linear constraints 
and we are consequently led to non-linear normal equations. It is 
seldom that such equations can be easily reduced to a linear set, and 
iterative solutions or other special methods are required. Similar 
difficulties arise with many hypotheses associated with problems of 
type D. 
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An alternative approach to questions of linearity of the change-over 
points might be via the fitting of lines to the observed unconstrained 
change-over points. As these will almost certainly be determined with 
different accuracies for different sets of data depending upon the numbers 
of points, and their positionings, in each phase, some weighting process is 
indicated, and the approach is entirely different from that adopted here. 

Finally, ‘from the biological viewpoint it must be emphasised that 
the ability to use techniques of the type suggested here in no way 
reduces the need for the experimenter to consider the biological reason- 
ableness of the model he uses, and whether it is more than a convenient 
approximation. 
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165 NOTE: A Practical Application of a 
Theoretically Inefficient Design’ 


C. P. Cox® 
National Institute for Research in Dairying 
University of Reading 
Reading, England 


1. INTRODUCTION 


In experiments on milking machine techniques for reducing the 
time taken to milk each cow, it was desired to compare milking times 
given by combinations of two pulsation ratios, P,; and P, with two levels, 
V, and V, , of vacuum in the machine line. It was expected (see Sec- 
tion 3) that the animal to animal variability would differ substantially 
for the different treatments, in view of which, and the short time avail- 
able for the experiment, it was decided to limit the experiment by 
having only two treatments per animal. The main restriction de- 
termining the design arose because the treatment combination P,V, 
was of such special interest that it was desired to obtain an estimate 
of the corresponding milking time from all the available animals. The 
simple solution used nicely illustrates the dual features of estimation 
and discrimination which often interplay in biological operational 
research. 


2. DESIGN PRINCIPLES 


Each animal formed an incomplete block receiving two of the four 
possible treatment combinations with four successive observations on 
each. The special feature is that, to provide the estimate required, 
all animals received the combination P,V, and practical considerations 
made it necessary for this to be the treatment for each animal in the 
first period. In theory therefore, the treatment differences would be 
confounded with any period differences but, in this case, previous 
experience had shown that these could safely be discounted in such a 
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short experiment. It was logical therefore to allot the three remaining 
treatment combinations equally, and randomly, to the animals in the 
second period. Thus, if P,V, is regarded simply as a control, the design 
itself is that properly deprecated by Yates [1936] for cases when treat- 
ment differences only are required, because of its inefficiency relative 
to the arrangement using all the six possible pairs. This latter design 
was inadmissible here because of the primary estimation requirement 
and it will appear that, because every observed difference contributes 
to each, the main effects of treatments and the interaction are estimated 
with sufficient accuracy to give discrimination at a useful practical level. 


3. ANALYSIS 


Let d,;; be the difference between the two period means for the jth 
animal on the treatment arrangement distinguished by the suffix 1, as 
follows; 


period 1 2 

t=1 PiV, — PV, 

1=2 — P,V2 

i=3 P\V, — 
and j = 1, --- , n for each 7. 


Let d;, be the mean of d;; over the n animals receiving the ith 
arrangement. Then the treatment comparisons are estimated as: 


d,,/2 d3./2 


main effect of P = 4(P2 — Pi)(Vi + V2) -1 1 -1 
main effect of V = 4(Pi + P:)(V2 — Vi) 1 -1 -1 


interaction PV = = Pi)(Va = Vi) 1 1 -1 


Hence, if o? is the variance of an individual d;;, each effect and the 
interaction is estimated with variance + + o3)/4n. 

With n animals in each of the three i-groups separate analyses of 
the d,; are obtained for each group according to the simple partition: 


treatment difference 1 
between animals (n — 1) 


to which a within animals term may be added if required. 
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In the example, the mean differences d, together with the between- 
animal mean squares with n = 8 were: 


treatment difference d;. mean square (7 d.f.) 
— —0.244 0.1467 
PiV, — PiV2 —0.638 1.5225 
Pi\V, — —1.631 1.9339 


We thus obtain the main effects and interaction, in minutes, as 
P = 3(0.244 — 0.638 + 1.631) = 0:619, 

and similarly, 
V=1.012, PV = 0.375. 


The appreciable differences between the above three mean squares 
are noteworthy and were expected because milking times are functions 
both of the machine variables and their interactions with highly in- 
dividual characteristics such as the amounts of milk secreted, neuro- 
hormonal responses and anatomical details of teats and udders. 

The variance of each main effect and of interaction PV is calculated 
as 


(0.1467 + 1.5225 + 1.9339) /32 


giving a standard error of 0.336 minutes. 

Since each mean square has the same number of degrees of freedom, 
the ordinary ¢-value for seven d.f. can be used for an approximate test 
of significance (after Cochran and Cox, [1957]) showing here that only 
the V-effect is significant at the five percent level. 


4. DISCRIMINATION VERSUS ESTIMATION 


It has been pointed out that the estimation requirement is fulfilled 
at the expense of the accuracy achieved for the treatment comparisons. 
Ilad a design employing all pairs of treatment combinations been 
applicable there would have been three more treatment groups, 


+= 4 —P,V., 
5 — P,V,, 
P,V,—P.V. 


The main effect of ?, for example, calculated without recovery of 
inter-animal information, is then, 
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and, instead of n, there will now be n/2 animals in each group. 

The variance heterogeneity in the particular example complicates 
the direct comparison but, if in a general case, transforming if necessary, 
we assume that each d;; has variance o”, the variance of a treatment 
comparison is now 


4 

16 n 
which compares with 307/4n obtained when only three differences are 
directly investigated. 

Both designs provide unbiassed estimates of the milking time on 
P,V, and the above loss of discrimination is compensated by the fact 
that this mean milking time is estimated with only half the variance 
which would have obtained if all the paired differences had been used. 

Finally, with reference to the particular example, when not only 
the mean but the distribution of milking times on the P,V, combination 
is of interest, the larger numbers of observations obtained with the 
design used will be additionally advantageous. 


= 0°/2n, 
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166 QUERY: On a Graphical Sequential Test 


Bross [1952] describes a graphical sequential test, and among the 
plotting rules is the instruction that “if new and old treatments lead 
to the same outcome then no information about superiority is obtained 
and nothing is plotted on the chart”. From his chart therefore dis- 
tinguishable results are considered significant whether they are 6 out 
of 6, or 6 out of 100 in the other 94 of which both treatments gave the 
same result. If one treatment kills 100 patients, and the other kills 
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94 and cures 6, normal probability ideas seem to suggest that there is 
only a 6 percent chance that there is some difference between the 
treatments, and not that it is less than 5 percent that there is no differ- 
ence. Is it wise to ignore some of the experimental results? Should 
the instruction be reworded “‘if the fact that the outcome of the new 
and old treatments is the same provides no information towards the 
particular problem being investigated, nothing is plotted on the chart: 
if it supports a hypothesis that the treatments are not different, then 
the experiment is plotted in the square adjoining the last plot diagonally 
away from the origin’’? 


REFERENCE 
Bross, I. D. J., [1952]. Sequential medical plans, Biometrics 8, 188. 


ANSWER: 


Although the query concerns a sequential procedure it applies 
equally well to the fixed sample procedure—the Sign Test—which is 
the parent of this sequential procedure so let us look at the simpler 
example. Both the Sign Test and its sequential analogue would be 
based on a study plan where one member of each patient pair is randomly 
allocated to treatment A or B (the second member receiving the re- 
maining treatment). When both patients in the pair have the same 
response (both live, both die), the “‘tie’’ is ignored insofar as the sig- 
nificance test is concerned. Thus in the example of the query the test 
would be based on the six pairs where the patient on treatment A 
survived and the patient on B dies (and the zero cases with the opposite 
results). The Mosteller-Corrected Sign Test would be (6—1)?/6=4.17. 
This exceeds 3.84 so the results would be significant at the 5 percent 
level. 

Is it wise to ignore some of the experimental results (eg. the “ties’’)? 
Let me make it clear that the “ties’’ would not be ignored in other 
phases of the analysis which deal with the clinical importance of the 
difference between treatments. The question, then, is whether “ties” 
should be ignored in dealing with the specific point: Is there any differ- 
ence in the responses to the two treatments? There are two technical 
approaches which will provide an answer to the question. First we 
might ask: Does ignoring the ties impair the efficiency of the procedure? 
The technical answer is that the loss of efficiency is negligible. Second 
we might ask: If significance tests are set up which involve the tied 
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observations, will they be better than the Sign Test? The technical 
answer is that there is no advantage to the alternative procedures and 
sometimes they result in a loss of power. While both technical ap- 
proaches lead to the answer ‘ties can be ignored”’ they may not satisfy 
an investigator who would like a simple rationale rather than a formal 
proof. 

Therefore I will indicate a rationale for the clinical situation. The 
purpose of therapeutic intervention is to change the natural clinical 
course of a disease in a favorable direction. Let us suppose treatment 
B fails to change the natural course (i.e. to death) while treatment A 
occasionally can change this course. Then the number of changes 
achieved on treatment A constitute the evidence of a difference between 
the two treatments. If there are 6 changes, we have the same amount 
of evidence of a difference whether it comes from 6 pairs or 600 pairs. 
If there were some spontaneous remissions (eg. not due to treatment), 
then we would have to compare the number of changes in the two 
series but this information’ would constitute the available evidence of a 
difference between the two treatments (irrespective of the number of 


pairs). When our interest lies in change of status then the “ties” are 
irrelevant. 


I. D. J. Bross 
Roswell Park Memorial Institute 
Buffalo, New York, U.S. A. 


167 NOTE: A Simple Method of Fitting the 
Regression Curve y = a + 6x + Bp” 


B. K. 
Department of Statistics 
University of Baroda, Baroda, India. 


Patterson [1956] described a simple method for fitting an asymptotic 
curve 


y=at Bp’, (1) 


for n = 4, 5, 6 and 7 equally spaced ordinates. In this method p is 
estimated by r, the ratio of two linear functions of y’s and a and 8 can 
be obtained by the linear regression of y on r’. 

This method can be extended to the problem of fitting the curve 
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y=at Bp" (2) 


for equally spaved ordinates. Here also the estimate of r is a ratio of 
two linear contrasts of the y's and the estimates a, 6 and 8 are calculated 
from a multiple linear regression of y on x and r’. 

As four parameters are to be estimated, the above method gives an 
exact fit for n = 4 using the formula 


Ya — 2y2 


j — 24, + Yo 8) 
The following formulae are suggested for n = 5, 6, 7 and 8: 


— — + Gyo ’ 


_ 10y, — 13y; — + 1d5y; 
10y5 13y.2 lly, + for n= 6, (5) 


— 10ys + ys — — 14ys 8y, + 


and 


10y; + — 5y, — 12y, — lly; — Ty2 + 20y, 
10y6 + — — 12ys — — Ty, + 


, forn = 8. (7) 


Then four formulae are chosen to give a reasonably high efficiency for 
independent y’s with equal variances. The percentage efficiencies for 
n = 5,6, 7 and 8, calculated from large sample formula for the variances 
of r and the efficient estimate (Shah and Patel, [1960]), are: 


p n=5 n=6 n=m7Z n=8 
0.0 95.2 91.3 88.4 86.4 
0.1 98.4 97.3 96.0 94.4 
0.2 99.7 99 .7 99.0 97.8 
0.3 99.9 99.7 98.4 96.9 
0.4 99.6 98.4 95.5 93.0 
0.5 99.0 96.1 91.4 87.5 
0.6 98.1 93.7 87.1 82.6 
0.7 97.3 91.4 83.0 76.1 
0.8 96.3 89.4 79.1 71.8 
0.9 95.4 86.7 75.9 68.2 


Thus the efliciencies of the proposed estimates are of about the same 
order as the efliciencies of Pattersons estimates of p for equation (1). 
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168 QUERY: On a Follow-up Study of the 


Growth of Children 


In conducting a long term follow-up study of the growth and de- 
velopment of a sample of children, the question arises as to whether 
losses between birth and 4 years were such as to bias the remaining 
- sample. Two specific questions may be asked: 
(a) Are losses independent of social class, age of parent etc? Leavers 
and stayers were compared, on several variables, using the Brandt- 
Snedecor method for a 2 X k table. 


Example 
Mother’s age Proportion 
at child’s Initial Stayers | remaining 
birth Sample | Leavers (X) (p) px 
up to 24 49 18 31 6327 19.6137 
25-32 117 25 92 7863 72.3396 
33+ 55 6 49 8909 43.6541 
Totals 221 49 172 p = .7783 | X& pX = 135.6074 
= 133.8676 
diff. = 1.7398 


x = (> pX — X)/pGg = 1.7398/.1725 
10.09 (.01 >P > .001). 


ll 


With several parental and social variables, x’ indicated significant 
differences between leavers and stayers. 

(b) What is the probability that the remaining sample (stayers) 
could arise by random sampling from the original sample? 

This would involve not a test of independence between staying and 
mothers’ age, but a test of the difference between observed frequencies 
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among stayers, and expected frequencies based on the proportions 
in the original sample. 

Can x’ be used to answer this question? It would seem illegitimate, 
in so far as expected values would themselves be derived from sample 
values, and not from population values. 

Douglas and Blomfield (1) employed a x’-test in this way, but in 
their case the initial sample consisted of all children born in a particular 
week in the majority of local authority areas of England and Wales. 
They derived expected frequencies in social class categories from the 
proportions in this sample, and applied them to the stayers. Is this 
legitimate, in so far as their expected frequencies might be said to be 
derived from population values? 

If the x*-technique is inapplicable in my case, is there any other 
way of assessing how unrepresentative the remaining sample has become? 


REFERENCE 


Douglas, J. W. B. and Blomfield, J. M. [1956]. The reliability of longitudinal sur- 
veys. Milbank Mem. Fund Quart. 34, 227-252. 


ANSWER: 


Since the original sample was random, the two questions, “Are 
the proportions of ‘leavers’ and ‘stayers’ effectively the same in the 
different rows?”’, and “Can the sample of stayers be regarded as a random 
sample of the original population?” are equivalent. The x’-test answers 
the question “‘Could the samples in all rows (or all columns) of the table 
have arisen by random sampling from a single population?’’, which is a 
slightly more specific version of the first question above. Since for 
your data x’ is significant, we are led to conclude that there is a genuine 
difference in the proportions of mothers of different ages in the sub- 
populations of ‘all stayers’ and ‘all leavers’ (i.e. of all mothers in the 
population who would have been classified under these headings had 
they been included in the sample): since the original sample was random, 
the sample of stayers must be regarded as having been drawn at random 
from the sub-population of all stayers, and so cannot be regarded as 
randomly drawn from the whole population. No independent test of 
this is necessary, or indeed possible: since the two questions are logicaily 
equivalent, the same test must be used to answer both. Admittedly 
it appears at first sight that a different test could be derived by com- 
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paring the observed numbers of stayers with their expectations, using 
the statistic 


2 
where X; is the number of stayers in the ith class, and the expectation 
and variance of X; are calculated on an appropriate null hypothesis. 
This statistic does not appear to depend on the numbers of leavers: 
however, since the true population proportions are here unknown, this 
test can only be based on the null hypothesis that the X; have arisen by 
independent random binomial sampling from the small numbers n, 
of the original sample, with a common probability estimated by 7. 
Hence Var (X;,) must be estimated by n,~@, so that formula (1) yields 


[X; np)” 


which, after a little manipulation, reduces to precisely the form that 
you have already used to calculate x’. Only if all the values of p are 
very small, so that g can be taken as effectively equal to 1, can formula 
(1) be approximated by the expression 


For your data g = 0.2217, so that the use of formula (2) would be 
illegitimate, as you conjecture, and would give a x’ less than one- 
quarter the size of the correctly calculated value. As far as I am able 
to judge from the information given in their paper, the test used by 
Douglas and Blomfield is in fact invalid for this reason. 

Having once decided that leavers and stayers do differ we may ask, 
quite legitimately, the further question, ‘(Is the non-randomness of the 
stayers, regarded as a sample of the original population, sufficiently 
serious to produce large biases in the estimates to be calculated from 
the sample?” If the overall proportion of leavers were small, the 
disturbances introduced by any variation in rate of loss might well be 
negligible in comparison with the original sampling variatien, so that 
the sample of stayers could be treated as though it were a random sample 
from the population: the resultant estimates would be biased, but the 
bias would be unimportant. If the true proportions in the various 
classes in the population are known, x’ may be calculated from formula 
(2): it must not be regarded as a test statistic, since the appropriate 
null hypothesis has already been rejected, but its size may be used as a 
criterion of the degree of non-randomness of the sample of stayers. 
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If the population proportions are unknown, as in your example, some 
indication of the seriousness of the risks involved in treating the sample 
of stayers as random may be obtained by considering the changes 
produced by leaving in the proportions of the sample falling in the 
various groups. For example, in the results quoted by Douglas and 
Blomfield the largest change in any employment group was in that of 
‘manual workers’, which comprised 0.691 of the original sample, and 
0.696 of the sample of stayers. Clearly, even though the authors 
were in error in maintaining that the sample of stayers could be regarded 
strictly as a random sample from the population, the degree of bias 
introduced by so regarding it was trivial. In your data, on the other 
hand, the proportion of mothers under 25 has fallen from 0.222 in the 
original sample to 0.180 in the sample of stayers, a much more serious 
change. 

Biases resulting merely from the changes in proportions can be 
largely avoided, in theory, by calculating separate estimates for the 
different groups, and combining them according to the proportions of the 
original sample. However, a breakdown of the sample taking into 
account all criteria that show significant changes in proportions may 
well leave you with groups too small for satisfactory estimation, and 
this method cannot allow for the further danger that the leavers will 
almost certainly be non-typical of the population even within groups, 
so that a considerable bias may remain even though the effects of 
disproportionality are successfully removed. 

M. R. SAmMPForD 
A. R. C. Unit of Statistics, 
Aberdeen, Scotland. 
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BOOK REVIEWS 


J. G. SKELLAM, Editor 
Members are Invited to Suggest Books for Listing or Review to the Editor 


16 DE JONGE, H. (editor). Quantitative Methods in Pharmacology. Amsterdam: 
North-Holland Publishing Company, 1961. Pp. xx + 391. £4.48. 


D. J. Finney, University of Aberdeen, Scotland. 


In May 1960, the Netherlands Region of the Biometric Society organized a 
Symposium for the Society. This volume contains the texts of the 26 papers. pre- 
sented, together with summaries of discussions and of remarks by chairmen. The 
papers differ greatly in interest and importance, and any selection for comment 
here may reflect the taste of the reviewer as much as the quality and balance of the 
Symposium. 

The first section of the book contains papers on sequential trials of drugs. 
HAJNAL presents an excellent account of sequential tests of four analgesics against 
a control, where assessment is in terms of a continuous variate and estimates of 
effects are wanted. JOHNSON introduces an important subject in a short paper on 
how to choose a sequential procedure that will meet practical requirements. RUMKE 
suggests uses of sequential techniques for estimating doses to be used in a definitive 
experiment; his idea is incompletely worked out, and questions of efficiency are not 
discussed, but it may give rise to improvements in the strategy of using pilot experi- 
ments. 

The second section, on drug standardization, includes a theoretical paper by 
CLARINGBOLD and EMMENS, consisting of a formal account of computational 
procedures appropriate to multivariate quantal responses under a general trans- 
formation, of which the normal equivalent deviate is a familiar example. Numerical 
illustration might have been helpful. Additional evidence is given that many types 
of data are insensitive to the choice of transformation, and the angle transform is 
therefore commended for general use. VAN STRIK’s suggestion of using rank as a 
metameter for semi-quantitative responses may sometimes be helpful, but the 
method fails completely in two special cases and one may suspect that other rank 
orders close to these extremes will give trouble. 

Evidently the most controversial session of the Symposium was that on non- 
parametric methods. HEMELRIJK and VAN DER VAART contribute papers 
on the efficiency and robustness of Wilcoxon’s test; both are of considerable theoreti- 
cal interest, though neither recognizes that the interpretation of numerical obser- 
vations can seldom be completed solely by appeal to tests of significance. VAN 
EEDEN’s ingenious but misleadingly titled paper illustrates this well; her method 
of estimating some points of the distribution function for a tolerance distribution, 
without assumption of any specified functional form, can evidently be extended to 
provide a method of estimating relative potency from two series of records which 
will work for some data, though for others it may be ambiguous. At best, it must 
be indeterminate to the extent of unit interval on the dose scale, and the author’s 
hopes of finding a method of assigning a measure of precision to such an estimate 
seems at present to lack foundation. JUVANCZ’s entertainingly vigorous attack on 
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non-parametric methods perhaps relies too much on what users of statistics have 
done rather than on what they ought to do, but comes closer to the needs of scientific 
research. In “‘A Visual Approach to Bio-Assay,’” GOLDBERG presents an ingenious 
system of graphs and scales that enable a close approximation to the maximum likeli- 
hood estimation of parameters for a quantal response regression to be achieved by a 
series of additions. The account of this, even for the single parameter exponential 
model, is so complicated as:to make this reviewer prefer the simplicity and speed of 
standard maximum likelihood calculations. 

The sections on drug screening and on mixtures of drugs are of consistently high 
standard. DUNNETT and SCHNEIDERMAN each discuss problems of evaluat- 
ing alternative screening rules with emphasis on economic aspects. Dunnett con- 
siders carefully the components of cost and loss that must enter into the specification 
of an optimal screening programme for a large number of drugs by appropriate 
decision rules; Schneiderman concerns himself more with the operating character- 
istics of his sequential procedures for screening anticancer drugs, and includes inter- 
esting comments on the use of a well-instructed computer. BROCK and SCHNEI 
DER seek to rehabilitate the therapeutic index by more careful definition than is 
customary; consideration of optimal experimental planning for the estimation of 
their index might throw more light on the seriousness of the loss in precision inherent 
in working at 5% and 95% response rates. Alternative, but not necessarily con- 
flicting, views on quantal responses to mixtures of drugs are discussed by ARIENS 
and SIMONIS and by HEWLETT and PLACKETT. The comprehensive mathe- 
matical models of the second will have the greater appeal for the statistician; for the 
pharmacologist, this summary of metrical aspects of Arién’s receptor theory should 
prove valuable, and perhaps the conjunction of the two will stimulate an effective 
synthesis. DIKSTEIN’s short account of a hypothesis relating chain length of 
primary and monoquaternary amines to excitation of guinea-pig muscle must surely 
be the forerunner of much work on the relation of physicochemical properties and 
pharmacological action. DE JONGH gives a useful and simple summary of the 
isobole method of representing different types of action of drugs in mixture. 

Inevitably a symposium of this kind includes a group of papers that defy any 
classification other than “miscellaneous’’, but those in the present volume are not 
the less interesting for that. GURLAND’s suggestions on bioassays for estimation 
of minute insecticidal residues include useful practical ideas for taking account of 
the masking of potency by inactive materials present in plant extracts and for the 
“fortification” of test extracts of very low potency so as to permit adequate assays 
when available quantities are severely limited. Finally, a contribution from 
BILLEWICZ describes a statistical study of a technique of measurement (the 
falling drop technique for estimating the concentration of deuterium oxide): al- 
though in itself purely physical, this is a good example of the manner in which a 
statistician can help to improve the accuracy of a laboratory technique. 


17 Le ROY, H. L. Statistische Methoden der Populationsgenetik. Stuttgart: 
Birkhauser Verlag, 1960, pp. 397, 67.5 Swiss Francs. 


C. Harte, Universitat zu Kéln, Germany. 


This book contains an extensive description of statistical methods used in 
modern population genetics. The reader is assumed to have some knowledge of 
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genetics and general biometrical methods, especially correlation methods and the 
analysis of variance. The first chapter deals with the influence of environment, 
genotype, and the combination of both on the expression of the phenotype. First, 
statistical and biometrical models are constructed. Then statistical methods 
are explained, and the use of the formulae demonstrated by numerical examples. 
The cases discussed are simple allelism, multiple allelism and polygenic inheritance. 
On this basis the next chapter gives the statistical treatment of the interdependence 
of genetical and environmental factors. Here the method of path coefficients is 
introduced, and by this method the analysis of correlations and interactions for 
different degrees of relationship is explained. : 

The third chapter deals with the same problems, but here they are analysed 
by means of the analysis of variance for the cases of a hierarchical and a two-way 
classification (incomplete and complete analysis of variance). The author gives both 
the method of analysis and the interpretation of the analytical results. The methods 
apply mainly to animal breeding, but in addition the results are generalised by dis- 
cussing the special problems which arise in plant breeding and human genetics. 

In the fourth chapter the methods given in Chapters 2 and 3 are applied to the 
special case of artificial selection in respect of one or more characteristics. Here not 
only the classical designs but also modern methods are considered. The author com- 
pares regression analysis, the method of path coefficients and the analysis of vari- 
ance. Different methods of treatment are also compared with respect to their effi- 
ciency and the amount of information they can give. 

In all chapters the formulae which are used are derived in a valid manner and 
the methods of analysis are illustrated by good examples. In several places there 
are diagrams which help to make the complex connections quite clear. The methods 
of analysis are clearly illustrated by detailed numerical examples. At the end of the 
book there is an extensive list of literature, which not only contains the papers 
mentioned in the text, but is in fact a bibliography on all the questions raised in the 
course of the discussions. In this way it is easy for the reader to find references in 
his own field of work on special questions which are not given a place in the book. 
The last pages give, besides the index and a list of authors, three very instructive 
tables. Their two-way classification according to problems and methods makes it 
possible to see the connections between the chapters and to locate the questions 
which are discussed under different topics from different points of view. 

The reader who’ knows enough in advance to follow the author and who has 
enough energy to work through nearlv 400 pages, will benefit greatly from this book. 


LI, C. C. Numbers from Experiments.—A basic Analysis of Variation. Pitts- 


os burgh: The Boxwood Press 1959. pp. 1x + 106, Cloth $4.00, paper $2.75. 


G. H. Jowett, University of Melbourne, Australia. 


This book is mainly concerned with explaining the algebra involved in the 
analysis of variance applied to one-way, hierarchical and two-way classifications, 
and to linear regression. 

The author is clearly fascinated by the partitioning of sums of squares and goes 
to considerable trouble to show how this happens. He makes considerable use of an 
interesting diagrammatic representation of a comparison between sums of squares 
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as a comparison between areas made up of ‘actual geometrical squares, which is 
apparently his own invention. 

Summation and suffix notations are used freely, particularly in later chapters, 
and matrices and orthogonal contrasts are also used, but only to demonstrate the 
splitting up of a sum of squares in an alternative manner. No basic statistical theory 
is given, nor is any serious attempt made to explain other aspects of the mathe- 
matical statistics. Some space is devoted to establishing expectations of mean 
squares, but this is based on quoted theorems of the type, 


+ y2 + + E(y2) + 
E(yy2) =u (yi, independent), 


the plausibility of which is justified mainly by reference to 1 very special case where 
y can take each value with equal probability. An attempt is also made to demon- 
strate the concept of a linear model; an ingenious example—‘Grandma’s birthday 
party”’—is used for this purpose, deviations from observed means being interpreted 
as gains and losses in a poker game by the grandchildren, the money being provided 
by Grandma with some interference by Grandpa. Illustrative examples used in the 
exposition are deliberately artificial, being chosen to employ simple numbers and 
thus concentrate attention on relationships and procedure rather than to suggest 
applications. 

Since the book does not cover the meaning and application of the techniques, 
it gives little help to anyone seeking to design or interpret the analysis of an experi- 
ment. On the other hand, it is not a suitable introduction to analysis of variance 
for a student of mathematical statistics or statistical method; it is loose and evasive 
about such important matters as degrees of freedom, testing hypotheses and inde- 
pendence, brings in sophisticated and inadequately explained jargon, and makes 
vague and unhelpful references to advanced mathematical topics (e.g. the rank of 
quadratic forms) in cases where the author has not been able to think of simple ex- 
planations. However, there is material here and there which, selected with dis- 
crimination, could give a teacher ideas for presentation or just tilt the balance of 
understanding for a student having difficulty with the algebra of the technique. 
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THE BIOMETRIC SOCIETY 


MEETINGS OF E, N. A. R. 


Annual Meeting, 1961 


The annual meeting of the Eastern North American Region of the Biometrics 
Society is being held during the meetings of the American Statistical Association 
at the Roosevelt Hotel in New York City, December 27-30, 1961. Jointly with 
the Biometrics Section of A. S. A., ten sessions of invited papers and one session 
of contributed papers have been arranged. 


Spring Meeting, 1962 

E. N. A. R. will meet jointly with the Institute of Mathematical Statistics at 
the University of North Carolina, Chapel, Hill, April 12-14, 1962. M. E. Turner, 
Division of Biometry, Medical College of Viginia, and A. F. Bartholomay, Biophysics 
Laboratory, Harvard Medical School, Cambridge, are co-chairmen of the E. N. A. R. 
program committee. H. B. Wells, Department of Biostatistics, University of 
North Carolina, Chapel Hill is handling local arrangements for E. N. A. R. 

Titles for contributed papers should be sent to M. E. Turner, who will then 
supply authors with abstract forms and appropriate instructions. ; 


Annual Meeting, 1962 


The annual meeting of E. N. A. R. in 1962 will be held during the meetings of 
the American Statistical Association at the Hotel Leamington in Minneapolis, 
September 7-10, 1962. Joint sessions with the Biometrics Section of A. S. A. are 
being arranged. 


W.N.A. R. 


MINUTES OF THE 1961 WNAR BUSINESS MEETING 


The annual meeting of WNAR, Biometrics Society was held at 2:30 P. M. 
June 16, 1961 at the University of Washington, Seattle. William F. Taylor, Presi- 
dent, WNAR, presided. 

The Secretary’s and Treasurer’s reports were read and approved. The WNAR 
student award (consisting of a free year’s subscription to Biometrics) was discussed 
at length and a motion was passed providing that any graduate student with a dual 
interest in the biological sciences and statistics should be considered eligible for the 
award. 

A possible meeting with the AAAS at Denver December 26-31 was discussed 
and the general concensus of opinion was that we should participate in this meeting. 

A motion was passed unanimously that WNAR meet for our annual meeting 
with AIBS at Oregon State University, Corvallis during the last week of August, 1962. 

The nominating committee gave their report and nominations were made from 
the floor. Those nominated were: for President, Carl A. Bennett, for the two 
Regional Committee vacancies, Gordon E. Dickerson, Charles Mode, Henry Tucker 
and Alvin Wiggins. 

The meeting adjourned at 3:20 P.M. Present at the meeting: Becker, Bennett, 
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G., Bennett, C., Berry, Bohidar, Brown, Calvin, Chapman, Hopkins, Leitch, 
Mantel, Nash, Petersen, Puri, Russell, Sandomire, Taylor, Tucker, Vaughan, 
Wiggins. 


WNAR ELECTION RESULTS 


Carl A. Bennett was elected WNAR President for 1962-63 and Charles J. Mode 
and Henry Tucker were elected to the Regional Committee for 1962-64. 


JOINT MEETING OF WNAR AND AAAS 


The Western North American Region will meet with the American Association 
for the Advancement of Science in Denver on December 28, 1961. There will be 
three sessions entitled ‘Experimental Design,” “Estimation of Populations (mobile 
and others),” and a ‘Contributed Papers Session.’”” Members wishing to present 
papers should contact Dr. Franklin A. Graybill, Statistical Laboratory, Colorado 
State University, Fort Collins, Colorado, Program Chairman. 


CHANGES OF MEMBERSHIP 
(July 15-October 15, 1961) 


Changes of Address 

Dr. Daniel J. Baer, 519 Chaucer Road, Dayton, Ohio, U.S.A. 

Mr. C. G. Barraclough, Chemistry Department, University of Melbourne, Parke- 
ville N. 2, Victoria, Australia 

Dr. Austin W. Berkeley, Department of Psychology, Boston University, 700 Com- 
monwealth, Boston 15, Massachusetts, U.S.A. 

Mr. Paul V. Blair, Institute of Enzyme Research, 1702 University Avenue, Univer- 
sity of Wisconsin, Madison 5, Wisconsin, U.S.A. 

Mr. Andre Bodeaux, Faculte Agronomique, Usumbura, Ruanda-Urundi, Belgian 
Congo 

Mr. Neeti R. Bohidar, Statistical Laboratory, Utah State University, Logan, Utah, 
U.S.A. 

Mr. Philippe F. Bourdeau, 370 Prospect Street, New Haven 11, Connecticut, U.S.A. 

Mr. N. Charliers, 199 Bonneville Sclayn, Province de Namur. Belgium 

Dr. J. B. Chassan, Hoffman-LaRoche, Nutley 10, New Jersey, U.S.A. 

Mr. Pan Leang Cheav, 1 rue Lucien Petit, Gembloux, Belgium 

Mr. William P. Chu, 5824 Fairfax Avenue, South Minneapolis, Minnesota. U.S.A. 

Dr. Ellsworth B. Cook, American Society Pharmacology and Experimental Thera- 
peutics, 9650 Wisconsin Avenue, NW, Washington 14, D. C., U.S.A. 

Mr. Tiberius Cunia, 2835 McCarthy Street, City of St. Laurent, Quebec, Canada 

Mr. R. N. Curnow, Unit of Biometry, The University, Reading, Berkshire, England 

Mr. Pierre T. A. Dagnelie, 53 Chaussee de Charleroi, Gembloux, Belgium 

Mr. Marc Dalebroux, 74 Planoy, Biesme-lez, Belgium 

Mr. R. de Coene, 79 rue Lincoln, Bruxelles 18, Belgium 

Prof. Martin E. Dehousse, 88 avenue Maeterlinck, Brussels 3, Belgium 

Dr. B. Diamantis, 2129 Scudder Street, St. Paul 8, Minnesota, U.S.A. 

Mr. H. M. Dicks, South Africa Sugar Association, Experiment Station, P. O. 
Mt. Edgecombe, Natal, South Africa 

Miss Martha W. Dicks, Division of Home Economics, Box 3354, University Station, 
Laramie, Wyoming, U.S.A. 


if 
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Mr. Victor A. Dirks, 12800 Dupont Avenue, S, Savage, Minnesota, U.S.A. 

Dr. N. R. Draper, Department of Statistics, University of Wisconsin, Madison, 
Wisconsin, U.S.A. 

Mr. Bruce A. Drew, The Pillsbury Company, 311 Second Street, S.E., Minne- 
apolis 14, Minnesota, U.S.A. 

Mr. L. Lee Eberhardt, Box 562, Concord, California, U.S.A. 

Mr. Harvey Eisenberg, 3722 Cedar Drive, Lochearn, Baltimore 7, Maryland, 
U.S.A. 

Dr. Khalil M. El-Kashlan, 16 Green Street, Moharrem Bey, Alexandria, Egypt 

Mr. Roger Firmin, 10 avenue R. Neybergh, Bruxelles, Belgium 

Mr. Andris Fogelmanis, 914 Sunset Drive, Pacific Grove, California, U.S.A. 

Mr. J. Fraselle, rue des Bas-Pres, Salzinnes (Namur), Belgium 

Dr. Seymour Geisser, Section on Biometry, Nat. Inst. of Arthritis and Metabolic 
Diseases, Bethesda, Maryland, U.S.A. 

Dr. Lue Goeminne, 6 Bagattenstraat, Gand, Belgium 

Mr. Richard A. Greenberg, 430 Central Avenue, New Haven, Connecticut, U.S.A. 

Mr. Samuel W. Greenhouse, Section on Statistics and Mathematics, National 
Institute of Mental Health, Bethesda, Maryland, U.S.A. 

Mrs. Fredrica G. Halligan, Warwick Towers, 9 Lafayette Court, Greenwich. Con- 
necticut, U.S.A. 

Dr. Emil L. Gumbel, Columbia University, 413 West 117 Street, New York 27, 
N. ¥., USA. 

Prof. John Gurland, U.S. Army Math. Research Center, University of Wisconsin, 
Madison. Wisconsin, U.S.A. 

Mr. M. A. Guzman, 10 Dixie Trail, Raleigh, North Carolina, U.S.A. 

Mr. A. M. Haider, Federation of Industries, Ottoman Bank Building, Baghdad, 
Iraq 

Dr. W. D. Hanson, Department of Genetics, North Carolina State College, Raleigh, 
North Carolina, U.S.A. 

Mr. Jacques Hardouin, 35 Via San Girolama, Perugia, Italy 

Mr. J. A. G. Hemptinne, 29 rue A. de Wasseige, Wepion, Fooz, Belgium 

Professor Dr. Leon Hennaux, Ville Beau Sejour, Coril Noirmont (Brabant) Belgium 

Mr. Jean Henry, 3 avenue Copyn Malaise, La Hulpe, Belgium 

Mr. John R. Howell, Route 5, Box 364, Orlando, Florida, U.S.A. 

Dr. Pierre O. Hubinont, 25 avenue Ad. Demeur, Bruxelles 6, Belgium 

Mr. Paul V. Hurt, 4721 Ames Street, Madison 5, Wisconsin, U.S.A. 

Mr. Floribert A. L. G. Jurion, 8 avenue Isidore Gerard, Bruxelles 16, Belgium 

Dr. Richard A. Lang, Department of Statistics, P. O. Box 5457, Raleigh, North 
Carolina, U.S.A. 

Dr. Lonnie L. Lasman, 124 Claridge, Eau Gallie 7, Florida, U.S.A. 

Dr. Jean Lebrun, 12 av. des Lucioles, Brussels, Belgium 

Mr. G. H. Lion, 61 rue Gabrielle, Bruxelles 18, Belgium 

Mr. Donald G. MacEachern, 1717 Brook Avenue, S.E., Minneapolis 14, Minnesota, 
U.S.A. 

Mr. John W. Mayne, 654 Sherbourne Road, Ottawa 3, Ontario, Canada. 

Mr. James H. Meade, Jr., Institute of Statistics, P. O. Box 5457, Raleigh, North 
Carolina, U.S.A. 

Mr. Hyman Menduke, 7808 Haines Road, Cheltenham, Pennsylvania, U.S.A. 

Mr. Forest L. Miller, 111 Villanova Road, Oak Ridge, Tennessee, U.S.A. 

Dr. Sigeiti Moriguti, Faculty of Engineering, University of Tokyo, Bunkyo-ku, 
Tokyo, Japan 
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Mr. Floyd R. Olive, McKamie, Arkansas. U.S.A. 

Mr. James G> Osborne. Old Chapel Hill Road, Route 1, Box 97, Durham, North 
Carolina, U.S.A. 

Mr. Jyun Otsuka, Kaname-machi 14, Toshimaku, Tokyo, Japan 

Miss Florence E. Petzel, Department of Home Economics, University of Texas, 
Austin, Texas, U.S.A. 

Mr. Andre Pieters, 39 Eeklastraat, Mariakerke (Gent) Belgium 

Mr. Jean Reuse, Clos des 4 Vents, Rue de Moorsel, (Brabant) Tervueren, Belgium 

Mr. Julien F. M. Ronchaine, 70 Grand Rue, Gembloux, Belgium 

Mr. Charles E. Rossiter, Pneumoconiosis Research Unit, Llandough Hospital, 
Penarth, Glamorgan, Wales 

Dr. Francesco Sella, Radiation Committee, Room 3468, United Nations, New 
York, N. Y., U.S.A. 

Mr. A. B. Siegelaub, 721 Walton Avenue, New York 21, N. Y., U.S.A. 

Mr. William A.-Small, P. O. Box 122-A, Tennessee Polytechnic Institute, Cooke- 
ville. Tennessee, U.S.A. 

Dr. H. Fairfield Smith, FAO-UN, P. O. Box 1555, Teheran, Iran 

Dr. John H. Smith, Department of Mathematics and Statistics, American Uni- 
versity, Washington 16, D.C., U.S.A. 

Mrs. Grace Scholz Spitz, 803 Marlo Drive, Falls Church, Virginia, U.S.A. 

Mr. Francois Sterckx, Aartcentrum 66, Geel, Belgium 

Dr. Howard L. Stier, United Fruit Company, 30 St. James Avenue, Boston 16, 
Massachusetts, U.S.A. 

Dr. Robert J. Taylor, CCNCS-NCI, National Institutes of Health, Bethesda, 
Maryland, U.S.A. 

Prof. Robert S. Temple, College of Agriculture, University of Tennessee, Knox- 
ville, Tennessee, U.S.A. 

Mr. George W. Thomson, 15093 Faust Avenue, Detroit 23, Michigan, U.S.A. 

Mr. M. Ulehla, 8 Carlysle Street, E. Hawthorn, Victoria, Australia 

Prof. A. van den Hende, 233 Coupure Links, Gand, Belgium 

Miss Cecile Van Kerchove, 13 rue Archimede, Brussels, Belgium 

Mrs. Pearl A. Van Natta, 1303 South Eudora, Denver 22, Colorado, U.S.A. 

Dr. B. Veen, van Tienhovenlaan 13, Velp, Netherlands 

Mr. Glen F. Vogel, 6600 Luzon Avenue, NW, Washington 12, D. C., U.S.A. 

Prof. Maurice Welsch, 99 Rue General Jacques, Embourg, Belgium 

Dr. R. Wette, M. D. Anderson Hospital, 6723 Bertner Avenue, Houston 25, Texas, 
U.S.A. 

Mr. Robert F. White, Science Information Lab., Smith, Kline and French, Phila- 
delphia, Pennsylvania, U.S.A. 

Dr. Manfred Woelke, Killicklaan 689, Les Marais, Pretoria, South Africa 


New Members 


At Large 


Mr. A. M. Haider, Federation of Industries, Ottoman Bank Building, Baghdad, 
Irag 

Dr. Helmut V. Muhsam, Professor of Statistics, Hebrew University, Jerusalem, 
Israel 


Australia 


Mr. P. J. Brockwell, Department of Statisties, University of Melbourne, Victoria, 
Australia 
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Mr. Peter D. Finch, Department of Statistics, University of Melbourne, Victoria, 
Australia 

Miss Betty Laby, Department of Statistics, University of Melbourne, Victoria, 
Australia 

Mr. J. S. Maritz, Department of Statistics, University of Melbourne, Victoria, 
Australia 

Miss Ilze Raudzins, 36 Boulder Avenue, Redcliffe, W. Australia 

Mr. K. V. Richardson, Department of Statistics, University of Melbourne, Victoria, 
Australia 


ENAR 


Mr. Charles Anello, 5905 Montgomery Street, Baltimore 7, Maryland, U.S.A. 

Mr. Charles E. Antle, 1620 University Avenue, Stillwater, Oklahoma, U.S.A. 

Mr. Thomas P. Bogyo, Department of Experimental Statistics, North Carolina 
State College, Raleigh, North Carolina, U.S.A. 

Mr. DuWayne C. Englert, Population Genetics Institute, Purdue University, 
Lafayette, Indiana, U.S.A. 

Prof. George A. Ferguson, Department of Psychology. McGill University, Montreal, 
Canada 

Mr. R. E. Foster. Forest Pathology Investigations, Department of Forestry, Vic- 
toria, B. C., Canada 

Mrs. Iris M. Kiem, Department of Medicine, University of Miami, Miami 36, 
Florida, U.S.A. 

Miss Joan Klebba, 1469 South 28th Street, Arlington 6, Virginia, U.S.A. 

Mr. James L. Pate, Department of Psychology, University of Alabama, University, 
Alabama, U.S.A. 

Dr. Norman B. Rushforth, Developmental Biology Center, Western Reserve 
University, Cleveland, Ohio, U.S.A. 

Mr. Jack Sebring, 1515 Ogden Street, NW, Washington 10, D. C., U.S.A. 


Germany 
Dr. O. Dittmar, Inst. f. Forstwiss, Alfred-Moeller-Str., Eberswalde, Germany 
Sweden 


Prof. Dr. E. Effenberger, (Hamburg 13) Gartenstrasse 2, Alsterchaussee, Erfurt, 
Sweden 


Switzerland 


Dr. Klaus G. Konig, Zahnarztliches Institut der Universitat, Zurich 28, Postfach, 
Switzerland 

Dr. Thomas M. Marthaler, Dental Institute, University of Zurich, Postfach, 
Zurich 28, Switzerland 

Dr. Heinrich Wagner, Ungargasse 32, Wien III, Austria 

Dr. K. Wuhrmann, Hofstr. 139, Zurich 44, Switzerland 


WNAR 


Dr. G. Eric Bradford, Department of Animal Husbandry, University of California, 
Davis, California, U.S.A. 

Mr. Howard F. Gingerich, 181714 Berkeley Way, Berkeley 3, California, U.S.A. 

Mr. Melville R. Klauber, Bldg. 315, Stanford Village, Stanford, California, U.S.A. 
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Dr. Orsell M. Meredith, Nuclear Medicine and Radiation Biology, UCLA Medical 
School, Los Angeles 24. California, U.S.A. 

Mr. William L. LeStourgeon, Statistics Group, Los Alamos Scientific Labs., Los 
Alamos, New Mexico, U.S.A. 

Dr. Daniel G. Siegel, Atomic Bomb Casualty Commission, Navy 3912, FPO, 

' San Francisco, California, U.S.A. 

Dr. Ervin P. Smith, Animal Science Department, Montana State College, Bozeman, 
Montana, U.S.A. 

Mr. N. Scott Urquhart, Computing Center CSURF, Colorado State University 
Fort Collins, Colorado, U.S.A. 
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NEWS AND ANNOUNCEMENTS 


Members are invited to transmit to their National or Regional Secretary (if mem- 
bers at large, to the General Secretary) news of appointments, distinctions, or retirements, 
and announcements of professional interest. 


NEWS ABOUT MEMBERS 


Jacob B. Chassan has been appointed statistician for the Research Department 
of Hoffmann-La Roche, Inc., at Nutley, N. J. 

C. Philip Cox, formerly of the N.I.R.D., Shinfield, England, has joined the staff 
of the Statistical Laboratory and Department of Statistics, Iowa State University, 
as Associate Professor. Professor Cox will be teaching and conducting research 
in the area of biological statistics. 

Seymour Geisser has been recently appointed Chief of Biometry of the National 
Institute of Arthritis and Metabolic Diseases. 

Dewey L. Harris, has been appointed as Assistant Professor to the staff of the 
Statistical Laboratory and Department of Statistics, Iowa State University. Dr. 
Harris will be consulting, teaching and doing research in the area of Genetics 
Statistics. 

Donald H. Hazelwood, formerly NIH Fellow, Department of Zoology, Washing- 
ton State University has taken a position as Assistant Professor of Zoology, Uni- 
versity of Missouri. 

Donald K. Hotchkiss, has been appointed as Assistant Professor to the staff 
of the Statistical Laboratory and Department of Statistics, Iowa State University. 
Dr. Hotchkiss will be consulting, teaching, and doing research in the area of statistics 
applied to animal nutrition. 

Robert Macey is now Assistant Professor, Department of Physiology, Uni- 
versity of California, Berkeley, moving from the University of Illinois School of 
Medicine. 

J. N. K. Rao of the Statistical Laboratory, and Department of Statistics, Iowa 
State University, has been promoted to rank of Assistant Professor. Dr. Rao wil! 
be teaching and conducting research in survey sampling in his new assignment. 

Elmer E. Remmenga, Chief, Computing Center, Colorado State University, 
will be on sabbatical leave from July, 1961 to July, 1962. 

Patrick K. Tomlinson has recently become an Aquatic Biologist, Ocean Salmon 
Investigation, California Department of Fish and Game. He formerly was ‘n the 
Biostatistics Section. 

Richard W. Vail, Jr. has taken a position as statistician in the Aerospace Cor- 
poration, Florida. He previously was statistician for Aerojet General Corporation, 
Azusa, California. 

Pearl A. Van Natta, formerly Biometrician at the Child Research Council, 
Denver, Colorado, is a self employed statistical consultant. 


SOUTHERN METHODIST UNIVERSITY—NEW PROGRAM IN 
STATISTICS 


Southern Methodist University announces the opening of a new Department 
of Mathematical and Experimental Statistics and Statistical Laboratory in the 
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Graduate School. Faculty members include Drs. Paul D. Minton and Vanamamalai 
Seshadri, assisted part-time by Mr. Del West and Mr. Mack Usher. 

Coursework leading to the degree of Master of Science will be offered, including 
courses in mathematical statistics, experimental statistics, probability, sampling 
theory and design of experiments. 

Departmental assistantships are available. 

For further information, write to Dr. Paul D. Minton, Department of Mathe- 
matical and Experimental Statistics, Southern Methodist University, Dallas 22, 
Texas. 


FLORIDA STATE UNIVERSITY—EXPANDED PROGRAM IN STATISTICS 


Beginning in September 1962, the Department of Statistics at the Florida 
State University will expand its graduate program to study and research leading 
to the Doctor of Philosophy degree in statistics. The success of the present program, 
currently limited to study to the Master of Science degree in statistics, has led to 
the approval of the new advanced training. The curriculum will be modified and 
expanded to include advanced work in statistical inference and decision theory, 
stochastic processes, advanced probability, multivariate analysis, operations re- 
search, special topics in biometry, theory of general linear hypotheses, non-parametric 
statistics and sequential analysis. Advanced courses will be offered in alternate 
years. The program leading to the Master of Science degree will be modified to 
permit both the thesis and non-thesis types of programs. 

The new program will be assisted through the recent addition to the staff of 
Dr. Vincent Hodgson from the London School of Economics and Political Science. 
Dr. Hodgson joints Drs. Ralph A. Bradley, Frank Wilcoxon, Richard G. Cornell 
and S. K. Katti in the Department of Statistics. Additional faculty appointments 
within the next two years are anticipated. 

Facilities have been greatly improved through completion of a new mathe- 
matical sciences building in which instructional, laboratory and office space was 
designed for use of the faculty and graduate students in statistics. The new building 
also houses an IBM 709 computer for research use. 

The Department has a limited number of teaching and research assistantships 
available. Proposals for three-year graduate fellowships are pending and such 
fellowships should be available for graduate students entering the new program in 
September. Interested students are invited to write to Dr. R. A. Bradley, Depart- 
ment of Statistics, the Florida State University, Tallahassee, Florida for further 
information. 
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CORRECTIONS 


Chin Long Chiang [1960]. A Stochastic Study of the Life Table and Its Appli- 
cations: I. Probability Distributions of the Biometric Functions. Biometrics 
16, 618-35. 


In formula (2), page 621, for 1. read I,. In formula (9), page 622, n should be 
deleted. 


Chin Long Chiang [1961]. A Stochastic Study of the Life Table and Its Appli- 
cations: III. The Follow-up Study with the Consideration of Competing Risks. 
Biometrics 17, 57-78. 


In line 1, page 62, m should be replaced by n. In formula (21), page 66, for 
S?_ read Sj, . In formula (32), page 69, for Qze 1 read Qskaz . In formulas (57) 
and (58), page 73, the super 2 should be deleted from o?. In formula (58), page 73, 
for Azz read Ax. In column (2) of table 4, page 77, for 835.87 read 935.87. In line 5 
from the bottom, page 77, for Qs read ex . 


J. K. Abraham [1960]. 154 Note: On an Alternative Method of Computing Tukey’s 
Statistic for the Latin Square Model. Biometrics 16, 686-691. 


In Table 3, cn‘ should multiply each }-‘, }>i and }>* appearing, and in Table 4, 
beneath “Treatments”, >-* should be similarly corrected. 


-— BOOKS OF IMPORTANCE IN YOUR FIELD 


Ready in January .. 
FIELD PLOT TECHNIQUE 


Completely revised, this edition is both a student text and a researcher's ref- 
erence manual. Practical aspects of statistical procedures, detailed examples 
of analysis are included . . . the book brings together and explains simply as 
possible the general principles, techniques, assumptions and guides to pro- 
cedures applicable to experimentation in agriculture and biology. 
Erwin L. LeClerg, United of 


Agriculture, 
rren H. Leonard and Andrew G ‘olorado State University, 
January 1962, about 350 pages, price pe 


An Introduction to the 
EXPERIMENTAL METHOD 


This book directs special attention to the design of experiments and the evalua- 
tion of variation by the less complicated statistical methods. The few technical 
terms which are used are fully explained. Mathematical manipulations are at a 
level of arithmetic and secondary school algebra, Both the reader who wants 
to understand research reports, and the investigator who will go on to more 
complex treatises will find the book a valuable aid. 


Y James Maxwell Little, Bowman Gray School of Medicine, Wake 
orest College, 1961, 88 pages, $3.00. 


order from 


BURGESS PUBLISHING COMPANY 


126 South Sixth Street ¢ Minneapolis 15, Minn. 
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Développement des méthodes biométriques et statistiques dans la recherche 
agronomique au Congo Belge et au Ruanda-Urundi.......... J. M. Henry 


Développement des méthodes biométriques et statistiques dans la recherche 
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variétés de betterave sucriére en Belgique........ M. Simon et N. Rousse. 


Les méthodes d’échantillonnage en sylviculture............. ANNE LENGER 


Radiosensibilité des graines d’andropogon issues de terrains uraniféres et non- 
uraniféres du Katanga...... Dr. J. MEewiIssEn, J. DaMBLON, ET Z. M. Baca 


Dix années d’activités de la Biometric Society en Belgique, au Congo Belge 
et au Ruanda-Urundi 


Revue bibliographique. 

Tome II, N° 1 Mars, 1961 

Quelques problémes de numérations bactériennes................ H. Grimm 

Comparaison de deux méthodes d’analyse de données non orthogonales dans 
le cas d’expériences 4 deux facteurs. ..... P. ET B. GRoSSMANN 

Analyse de données non orthogonales dans le cas d’une expériences 4 deux 


Revue bibliographique: D. J. Finney: Experimental Design and its Statistical 
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JOURNAL OF THE 
AMERICAN STATISTICAL 
ASSOCIATION 


Volume 56 December 1961 Number 296 


TABLE OF CONTENTS 


Mectitying Tnepection: OF F. J. ANSCOMBE 
On Comparing Intensities of Association between two Binary 
Characteristics in two Different Populations ............. Aanes 


Occupational Components of Educational Differences in Income 
Otis DupLey Duncan 
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1757 K Street, N.W., Washington 6, D. C. 


7 
318 
Bis 
- i} 


TECHNOMETRICS 
A Journal of Statistics for the 
Physical, Chemical and Engineering Sciences 


Vol. 3, No. 3 August, 1961 


CONTENTS 
The 2"? Fractional Factorial Designs....... G. E. P. Box anp J. S. HunTER 


Reduced Designs of Resolution Five....J.C. WHITWELL anp G. K. MorsBey 


Partial Confounding in Fractional Replication.............. W. J. YoupEen 
Finding New Fractions of Factorial Experimental Designs...... R. E. Fry 
A Study of the Group Screening Method.................... G.S. Watson 
Missing Values in Response Surface Designs........... Norman R. Draper 


Use of Tables of Percentage Points of Range and Studentized Range 
H. Leon Harter 


The Optimum Allocation of Spare Components in Systems 
Donatp F. Morrison 


The Reliability of Components Exhibiting Cumulative 


An Analysis of Some Relay Failure Data from a Composite 
Exponential R. R. Prairie anp B. Ostie 


Applications of Truncated Distributions in Process Startups 
H. SmituH anp D. W. Grace 


Estimating the Poisson Parameter from Samples that 
Are Truncated on the A. Cuirrorp Conan, Jr. 


A General Simulation Programme for Batch Chemical Plants 
J. Dyson, P. L. GotpsmitH anp J.S. M. Rosertson 


Statistical Programs for High Speed Computers 


Notices 


Technometrics is published quarterly in February, May, August and 
November. To members of the American Statistical Association and the 
American Society for Quality Control the annual subscription rate is $6.00. 
The annual non-member subscription rate is $8.00. Checks should be made 
payable to Technometrics and addressed to Technometrics, Post Office Box 
587, Benjamin Franklin Station, Washington 6, D. C. 


| 
SS 
| 
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The aim of this journal of abstracts is to give complete coverage 
of published papers in the field of statistical theory (including 
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and new published contributions to statistical method. 


All contributions in the following five journals—being wholly 
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organolepsis 
Successive approximation, see iteration 
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Sufficient statistic, 2, 121, 585, and see 
efficiency, estimation 
Survival curve, 173, and see dose 
response curve, mortality, time 
response curve 
Survival time, see time response curve 
Symmetrical complementation design, 
499 
Synergism, 188 
Tables, miscellaneous, 52, 55, 91, 92, 
100, 126, 244, 245, 247, 291, 321, 
324, 408, 495, 496, 501, 507, 533, 
534, 535, 544, 546, 561, 593, 600, 652 
Taste test, see organolepsis 
Taxonomy, 331 
Taylor series, 122 
Test, see analysis of variance, chi 
square, comparison, confidence 
limits critical region, goodness of 
fit, order, rank, rejection of data, 
sequential, sign test, triangle test, 
t test 
asymptotic, 142 
distribution-free, 182 
exact, 133 
conditional, 595 
multiple comparison, 176, 324, 539 
multiple range, 
Duncan’s, 321, 540 
Newman-Keuls, 540 
Tukey’s, 540 
multivariate, 496, 541 
nonparametric, 183, 539, 657 
of boundedness, 166 
of distribution, 507 
of significance, 
of difference between regressions, 
637 
largest difference, see multiple range 
test 
power of, 176, 501 
sequential, 176, 178, 179, 622 
graphical, 649 
using range, see multiple range test 
Wilcoxon’s, 543 
power of, 182 
robustness of, 183 
Theory, see biometry, hypothesis, law, 
model 
Therapeutic index, 186, 658 


Tie, 649 
Time response curve, 634, and see dose 
response curve 
Time series, 391 
Toxicology, 186, 188, 501, 658, and see 
bioassay, pharmacology 
Tracer, 189, and see radioactivity 
Transformation, 1, 180, 209, and see 
additivity, analysis of variance, 
bioassay, discriminant function. 
logistic curve, matrix, model 
angular, 5, 372- 
arc sine, see angular 
inverse, 1 
Laplace, 350 
logarithmic, 554 
logit, 1, 372, 501 
log log, 1 
power, 498 
probit, 94, 372, 501 
rankit, 251 
Trend, see regression 
Triangle test, 251 
Truncation, 57 
T test, 181, and see confidence limits, 
Hotelling’s T 
multiple, 324 
non-normal data, 182 
power of, 182 
Tuberculin, 33, and see Mantoux 
Tumor, see cancer 
Variance, see covariance, genetic 
adjusted, 112 
analysis of, 258, 450, 562, 659, and 
see additivity, analysis of co- 
variance, analysis of groups of 
experiments, chi square, compo- 
nents of variance, error, fitting 
constants, interaction, interblock 
error, least squares, model, multi- 
variate analysis, path coefficients, 
regression, structural analysis, 
test, transformation 
model II, 175 
of differences, 164 
random effects, 608 
approximate, 95, 368 
asymptotic, 58, 558 
components, 145, 438, 490, and see 
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components of covariance, cumu- 
lant components, _ structural 
analysis 
computation of, 607 
genetic, 233 
variance of, 361 
error, 506 
estimated, 273 
genotypic, 232 
homogeneity of, 90 
interblock, see information 
matrix, see covariance matrix 
of difference of adjusted means, 116 
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of difference of means, 452 
of genetic correlation, 359 
of product, 172 
of proportion, 31 
of regression coefficient, 286 
of scale parameter, 503 
proportional to power of mean, 498 
Vital statistics, 167 
Weibull, see distribution 
Weighting, 2, 122, 300, 481, 581 
choice of, 288 
Working value, 2 
X-ray, see radiology 
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BACK ISSUES 
Back issues of Biometrics are available at the following postage-paid 
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